ADA131275 


ONA  6153F 


ALTERNATE  IMPLOSION  MODELS  FOR  THE 
PLANE  PARALLEL  DIODE 


Robert  E.  Terry 
John  U.  Guillory 
JAYCOR 

205  South  Whiting  Street 
Alexandria,  Virginia  22304 

2  April  1982 


Final  Report  for  Period  1  November  1980—1  November  1981 


CONTRACT  No.  ONA  001-79-C-0189 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED. 


DTIC 

#%ELECTE 
W  AUG  1  1  1983 

HiB 


THIS  WORK  WAS  SPONSORED  BY  THE  DEFENSE  NUCLEAR  AGENCY 
UNDER  RDT&E  RMSS  CODE  B323082466  T99QAXLA00020  H2590D. 


Prepared  for 
Director 

DEFENSE  NUCLEAR  AGENCY 
Washington,  DC  20305 


88  08  11  008 


Destroy  this  report  when  it  is  no  longer 
needed.  Do  not  return  to  sender. 

PLEASE  NOTIFY  THE  DEFENSE  NUCLEAR  AGENCY, 
ATTN:  STTI ,  WASHINGTON,  D.C.  20305,  IF 
YOUR  ADDRESS  IS  INCORRECT,  IF  YOU  WISH  TO 
BE  DELETED  FROM  THE  DISTRIBUTION  L*3T,  OR 
IF  THE  ADDRESSEE  IS  NO  LONGER  EMPLOYED  BY 
YOUR  ORGANIZATION. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  This  PAGE  (Whan  Dm  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  2.  GO VT  ACC ESS!ON  NO^ 

DNA  6153F  /J/3|27to 

^RECIPIENT'S  CATALOG  NUMBER 

4.  TITLE  (mtd  Sutfltli) 

ALTERNATE  IMPLOSION  MODELS  FOR  THE 

PLANE  PARALLEL  DIODE 

S.  TYPE  OF  REPORT  •  PERIOD  COVERED 

Final  Report  for  Period 

1  Nov  80—1  Nov  81 

S.  PERFORMING  ORG.  REPORT  NUMBER 

J207-82-002 

7.  authorc*; 

Robert  E.  Terry 

John  U.  Guillory 

S.  CONTRACT  OR  GRANT  HUMBERT*.) 

DNA  001-79-C-0189 

S.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

JAYCOR 

205  South  Whiting  Street 

Alexandria,  Virqinia  22304 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Subtask  T99QAXLA000-20 

It.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

Director 

Defense  Nuclear  Agency 

Washington,  DC  20305 

12.  REPORT  DATE 

2  April  1982 

13.  NUMBER  OF  PAGES 

90 

MT  MONltoRING  AGENCY  NAME  *  AOORESSfU  dtlt*rm>t  tram  Controlttn*  OlHc*) 

IS.  SECURITY  CLASS  (ot  thi*  report; 

UNCLASSIFIED 

Wa.  DECL  ASSI  FlCATION/ DOWN  GRADING 

N/jfsTnce  UNCLASSIFIED 

ft.  DISTRIBUTION  STATEMENT  (o>t  Hil«  Report) 

Approved  for  public  release;  distribution  unlimited. 

17.  DISTRIBUTION  STATEMENT  (at  the  abelract  inNfM  In  Black  30,  It  dlflmnt  tram  Report) 

is.  supplementary  notes 

This  work  was  sponsored  by  the  Defense  Nuclear  Agency  under  RDT&E  RMSS 

Code  B323082466  T99QAXLA00020  H2590D. 

IS.  KEY  WORDS  (Centlnm*  an  rarer**  aide  It  neceeeery  and  Identity  Or  block  mmb*r) 

Plasma  Physics  Hydrodynamic  Codes 

Imploding  Plasma  Loads  Modeling 

Radiative  Sources  Implosion  Model 

Photon  Sources  Lagrangian  Model 

SO.  ABSTRACT  (Can tint*  «  rararc*  *10*  It  naca**arr  mtd  Identity  br  black  number) 

.The  evolution  of  imploding  Z-pInch  plasmas  in  a  p 

lane  parallel  diode  is 

PR,  described  here.  The 
d  solution  via  the  appropri- 
tion  coupled  hydrodynamic 
second,  in  an  approximate 
field,  combined  with  the 
ynamic  solution  is  „ 

^modeled  in  two  ways  through  the  computer  code  ZDI 
first  consists  In  a  complete  electromagnetic  fiel 
ate  Hertz  vector  potential,  combined  with  a  radia 
solution  of  the  Braglnskii  fluid  equations.  The 
electrodlffuslve  solution  for  the  axial  electric 
radiation  coupled  hydrodynamic  model.  The  hydrod 

DO  I  '22*73  1473  EOlTION  OF  1  NOV  SS  IS  OBSOLETE  UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  FACE  (Whm  Dm*  tnian« 


20.  ABSTRACT  (Continued) 


^  effected  by  an  implicit  advance  of  ion  temperature,  electron  temperature 
plus  chemical  potential,  and  (in  the  electrodiffusive  limit)  the  axial 
electric  field.  This  Is  subcycled  under  the  control  of  either  a 
predictor/corrector  or  fully  implicit  advance  of  the  radial  fluid  velocity 
and  Lagranglan  position  coordinate.  Radiation  transport  is  calculated  by 
probability  of  escape  algorithms  operating  on  a  preselected  set  of  frequency 
domains.  The  electromagnetic  field  advance  is  effected  through  wave  equation 
quadrature  (for  Z  and(31'Z)  followed  by  a  numerical  differentiation  of  this 
potential  (smoothly  Interpolated  to  the  fluid  mesh)  which  provides  and  B0. 


^  b  b  ' 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  of  This  FAGBflWiMi  Dm*  Snt*r*ri; 


TABLE  OF  CONTENTS 


Section  ?a§e 

PREFACE . ni 

I.  INTRODUCTION  .  1 

A.  Motivation .  1 

B.  Computational  Approach  .  6 

II.  ACCURATE  METHODS  FOR  NUMERICAL  DIFFERENTIATION  .  13 

A.  Survey  of  the  Theory  of  Smooth  Interpol  ants . 13 

B.  A  Family  of  Generalized  Splines  R(x,y) . 17 

C.  The  Interactive  Interpolation  Package 

(VAX  Implementation) . 20 

D.  Applications  of  Interest  in  Modeling  Diode 

Impl  oded  Pi  asmas . 22 

III.  THE  HERTZ  DRIVEN  IMPLOOING  PLASMA  RADIATOR 

(PROGRAM  ZDIPR)  .  27 

A.  Overview . 27 

B.  Theoretical  Sumnary . 23 

C.  Electromagnetic  Algorithms  and  Subroutines  .  32 

D.  Transport  Coefficients  .  34 

E.  Fluid  Density  and  Flow  Field  Time  Derivatives.  ...  37 

F.  Fluid  Temperature  Time  Derivatives  .  42 

G.  The  Subcycle  Calculating  Thermal  and  Electro- 

diffusive  Evolution  in  the  Model  Plasma . 48 

H.  The  Complete  Fluid  Evolution  Package  -  HYDROPUSH  .  .  55 

I.  The  Planned  Sequential  Benchmark  of  This  Method.  .  .  65 


1 


TABLE  OF  CONTENTS  (Continued) 


Page 


Appendix  I.  VAX  Interactive  Interpolation  Package . 59 

Appendix  II.  Derivation  of  the  Electrodlffusive  Limit  .  75 

Appendix  III.  Area  Weighted  Differences . 79 


PREFACE 


We  wish  to  gratefully  acknowledge  the  Important  and  central  contributions 
to  this  effort  by  J.  Davis,  P.  C.  Kepple  and  D.  Duston  at  NRL  in  providing 
fundamental  equation  of  state  and  radiation  emission  data  as  well  as  timely 
and  cogent  advice  throughout  the  developmental  phase.  We  are  also  grateful 
for  the  painstaking  and  effective  aid  given  by  J.  Apruzese  and  R.  Clark  at 
NRL  in  understanding  and  using  the  radiation  transport  algorithms  discussed 
here. 

In  addition  the  program  development  help  of  D.  W.  Mack,  S.  SI  inker  and 
L.  Roelofs  at  JAYCOR  Is  gratefully  acknowledged,  as  are  the  several  sub¬ 
routines  of  the  IMSL  library  used  here. 


I.  INTRODUCTION 


A.  Motivation 

1  2 

Existing  models  of  imploding  wire  or  puff  gas  plasmas  ’  leave  a 
number  of  physical  processes  unexplored  or  only  partially  examined  in  their 
assessment  of  diode  Imploded  plasmas  as  radiation  sources.  The  one-dimension¬ 
al  Implosion  code  to  be  described  here  is  an  attempt  to  increase  the  physical 
relevance  of  such  models. and  to  enhance  the  versatility  and  robustness  of  the 
nunerical  schemes  used  to  implement  them. 

Among  the  relatively  unexplored  physics  questions,  one  finds  those 
of  photon/plasma  energy  exchange,  plasma  thermoelectric  radial  stresses  and 
the  effect  of  separately  evolved  ion  and  electron  temperatures  on  opacity 
calculations.  The  incompletely  examined  physics  involves:  the  processes  of 
current  penetration  into  the  plasma  load,  the  correlations  of  diode  waveforms 
and  plasma  motions  (generated  by  reflected  power  in  the  diode  cavity),  the 
energy  partition  between  thermal  motion  and  ionization  potential,  and  the  role 
of  marginally  stable  drift-speed  limitations  in  determing  the  implosion  dynamics. 
These  considerations  can  be  organized  into  several  concrete  goals  for  the 
calculation. 

1)  Find  how  important  are  the  details  of  photon/plasma  energy 

exchange  in  structuring  the  implosion  trajectory  and  (spatially 
resolved)  energetics.  The  transitions  in  qualitative  characteris¬ 
tics  (ranging  from  refri gerati ve  collapse  to  strong  reflection) 
noted  in  the  core-corona  model"*  were  quite  sensitive  to  the  opacity 
calculation.  In  a  full  1-0  implosion  calculation  the  resolution 
of  the  enhanced  "thermal  conduction"  due  to  direct  photon/plasma 
energy  exchange  will  provide  a  more  accurate  picture  of  the  plasma 
•mission  profile. 
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2)  Determine  the  role  of  separate  ion  and  electron  temperatures 
in  establishing  this  emission  profile  and  in  governing  the  plasma 
dynamics.  As  compared  with  a  single-temperature  picture,  one 
expects  the  proper  two-temperature  treatment  to  alter  the  pressure 
gradients  in  low  density  regions  and  change  the  opacity  in  higher 
density  regions,  perhaps  altering  the  spectral  mix  of  line  and 
continuum  radiation. 

3)  Examine  the  effect  of  the  chemical  potential  (properly  included 
in  the  plasma  energy  balance  as  a  function  of  electron  energy  and 

4 

Ion  density)  on  the  Implosion  energetics  and  the  emission  profile. 

An  accurate  partition  of  Incoming  energy  will  probably  provide  lower 
(more  physical)  temperatures  at  peak  implosion  than  earlier  models. 

4)  Develop  an  accurate  picture  of  the  current  penetration  mechan¬ 
isms.  Diffusive  transport  of  Ez(r,  t)  (or  E^)  is  dominant  to  a  first 
approximation,  but  overall  Implosion  performance  is  sensitive  to  the 

e 

effective  skin  depth.  The  role  of  the  free  wave  field  components 
in  setting  this  skin  depth  has  not  been  examined,  and  the  rate  of 
Inward  propagation  for  the  drift- speed  limitations  in  Jz  may  be 
affected  by  these  components. 

5)  Resolve  the  sources  of  reflected  power  within  the  plasma  load 
and  relate  them  to  the  available  measured  waveforms  from  experiment. 

An  electromagnetic  calculation  of  all  fields  within  the  plane  parallel 
diode  opens  one  more  channel  of  diagnostic  information,  and  may 
provide  a  means  of  inferring  plasma  motion  from  measured  waveforms. 

It  will  also  provide  more  ac ’  irate  statement  of  energy  conservation 
because  the  energy  con*-  neo  in  the  net  free  wave  component  will  be 
Included. 
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6) 


Examine  the  effect  of  radial  ambipolar  electric  fields,  and  the 
resulting  stress,  on  the  implosion  dynamics,  particularly  in  the  low 
density  regions  where  drift  speed  limits  are  expected.  The  radial 
electric  fields  arise  from  thermoelectric  effects  and  axial  drifts 
of  charge  carriers.  Large  temperature  gradients,  allowed  by  weak 
cross  field  thermal  conduction  and  enhanced  ohmic  heating  in  the 
corona  plasma,  may  provide  a  significant  radial  electric  field  and 
resulting  Inward  stress  on  the  plasma. 

Once  these  physical  questions  are  addressed  coherently,  a  smooth 
upgrading  in  the  sophistication  of  the  radiation  emission  and  transport 
model  becomes  a  fruitful  exercise.  With  the  radiative,  hydrodynamic,  and 
electrodynamic  models  structured  from  first  principles  there  are  no  "free 
knobs"  and  the  subsequent  benchmarks  with  an  appropriately  configured 
experiment,  if  successful,  provide  a  firm  basis  for  similar  theory  in  more 
difficult  and  perhaps  more  practical  diode  configurations. 

Apart  from  these  physical  questions  are  various  topics  arising  in 
the  numerical  implementation  of  whatever  implosion  model  is  selected.  First 
is  the  choice  of  grid;  a  typical  simulation  will  show  the  model  plasma 
compressed  several  orders  of  magnitude  during  the  course  of  a  calculation. 

A  strict  Eulerian  fluid  calculation  can  thus  lose  its  spatial  resolution 
at  the  final  collapse,  when  this  resolution  is  most  needed.  Moreover  large 
density  gradients  are  commonplace,  and  it  is  probable  that  the  dynamics  of  the 
low  density  region  is  dominated  by  the  balance  of  density  gradients  and 
magnetostrictive  stresses.  A  Lagrangian  description  is  thus  preferable  in 
an  implosion/explosion  calculation,  so  long  as  the  external  stresses  on  the 
system  can  be  defined  readily  over  the  spatial  domain  of  interest.  On  the 
other  hand  the  fundamental  field  ( Ez }  generating  the  implosion  is  defined 
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in  the  fixed  laboratory  frame  and  is  most  naturally  discretized  on  an  Eulerian 
mesdi.  The  electromagnetic  problem  therefore  requires  an  interpolation  capa¬ 
bility  which  can  map  external  stresses  onto  a  convecting  fluid  and  extract 
field  sources  from  this  moving  fluid.  The  natural  compromise  adopted  here 
is  a  completely  Lagrangian  fluid  model  convecting  through  a  laboratory- 
frame  (Eulerian)  electromagnetic  grid.  The  electromagnetic  grid  can  be 
made  sufficiently  non-uniform  to  concentrate  the  Ez  information  density  near 
the  axis  while  the  smoothly  distorted  Lagrangian  fluid  grid  is  a  natural 
choice  for  an  accurate  resolution  of  density  gradients. 

A  second  concern  Is  the  choice  of  Ez  or  Be( Jz )  as  the  fundamental  electro 
dynamic  variable.  In  our  view  the  proper  selection  is  Ez  because  it  is  the 
natural  field  from  which  to  establish  drift-speed  limitations  on  the  current 
density.  This  consideration  plays  a  fundamental  role  in  both  the  development 
of  the  Hertz  vector  formalism  and  in  the  choice  of  the  diffusive  limit 
discussed  below  (III. A.).  Using  8g  and  0Z  *  *  Be)z  a  drift  speed 

limitation  algorithm  has  two  disadvantages.  The  current  is  derived  from  a 
potentially  noisy  differentiation  process.  The  limit  Jz  *  enfi  C$  does  not 
imply  a  completely  local  value  for  Ez  because  the  conductivity  depends  on 
Ba,  and  hence  on  at  other  grid  points.  It  is  clearly  possible  to  remove 
the  second  problem  by  Iteration,  but  the  physical  process  providing  the  drift 
speed  limit  Is  a  local  (fine  scale)  one.  By  proceeding  first  from  the  local 
values  of  Ez»  then  iteratively  establishing  small  corrections  to  the  (possibly 
limited)  Jz  from  non-local  considerations,  one  achieves  a  more  natural  and 
probably  less  noisy  convergence  to  the  limit. 

A  third  point  is  the  selection  of  forward  time  integration  methods. 

A  common  preference  for  conservatively  advanced  explicit  integrators  must 
be  examined  carefully  in  the  context  of  this  radiatively  coupled  problem. 


As  the  level  of  sophistication  in  the  radiative  transfer  physics  is 
Increased  a  considerable  expense  is  incurred  for  each  fluid  time  deriva¬ 
tive.  This  expense  will  scale  linearly  with  the  number  of  frequency  domains 
and  quadratically  with  the  spatial  resolution.  In  order  to  optimize  the 
complete  algorithm- it  is  therefore  necessary  to  make  the  best  possible  use 
of  every  time  derivative  available  from  the  space-time  mesh.  Moreover,  as 
distinct  from  calculations  dominated  by  processes  internal  to  the  fluid,  the 
problem  of  simulating  an  electromagnetically  driven  implosion  requires  a 
self-consistent  time  evolution  of  both  field  and  fluid.  This  class  of 
problems  is  more  likely  to  be  handled  efficiently  and  accurately  by  an 
implicit  method,  or  at  least  by  a  predi ctor/ corrector  scheme.  This  choice 
has  only  the  disadvantage  of  being  more  cumbersome  to  implement,  and,  as 
described  below,  the  methods  chosen  make  a  smooth  transition  from  "nearly 
explicit"  to  "fully  implicit."  As  dictated  by  the  performance  characteris¬ 
tics  of  the  complete  algorithm,  the  user  will  be  able  to  optimize  the 
accuracy  and  cost  tradeoff  by  several  means. 

A  final  consideration  is  the  accuracy  of  spatial  differentiation 
within  either  the  hydrodynamic  or  electrodynamic  algorithms.  If  these 
algorithms  can  make  good  use  of  non-uniform  meshes  in  space  and  time,  then 
one  can  sustain  useful  accuracy  with  fewer  points.  The  coimion  differentiation 
methods  based  on  finite  differences  are  quite  accurate  when  a  dense  mesh  is 
admissible,  but  if  one  wishes  to  keep  down  the  cost  of  radiative  transfer 
calculations  a  higher  order,  smoother  scheme  is  preferable.  Some  useful 
alternatives  are  discussed  in  Chapter  II. 
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B.  Computational  Approach 


The  specific  physical  content  of  the  model  has  been  described 
6  7  3 

elsewhere.  *  ’  The  task  of  implementation  can  be  readily  broken  down  into 
the  three  problem  areas  treated  below  -  electromagnetics,  hydrodynami cs , 
and  radiative  transfer.  The  requirements  for  each  problem  area  and  a 
summary  of  the  algorithms  employed  are  given. 

The  electrodynamic  problem  statement  is  simple:  given  voltage  wave¬ 
forms  at  the  periphery  of  the  diode,  establish  the  E  and  B  fields  on  the 
Interior  spatial  domain  for  the  duration  of  the  implosion  process.  By 
driving  the  diode  E/M  model  at  one  radius  with  a  voltage  waveform,  the  current 
drawn  by  the  load  and  the  voltage  waveforms  at  any  interior  point  are  avail¬ 
able  as  predicted  observables/diagnostics  of  the  code  and  the  experiment. 
Conventional  approaches  replace  the  vacuum  portion  of  the  diode  with  a 
circuit  relation  derived  from  Faraday's  law,  but  a  complete  electromagnetic 
calculation  provides  more  information.  If  this  later  goal  is  chosen  then  the 
Issue  of  efficient  and  effective  calculation  becomes  an  important  question. 
With  the  implosion  model  discussed  here  the  choice  has  been  to  admit  either 
option,  viz.  the  full  electromagnetic  computation  or  the  electrodiffusive 
limit  and  circuit  equation. 

For  the  first  option,  it  is  instructive  and  in  many  ways  practical 
to  utilize  the  generalized  Hertz  vector  potential.  First,  this  new  potential 
reduces  the  electrodynamics  to  a  single  component  of  a  single  vector  wave 
equation.  This  reduction  simplifies  the  calculation  by  placing  no  self- 
consistency  constraints  on  the  time  evolution  of  the  fields  E„,  B,.  The 
required  self-consistency  is  automatic  when  Ez,  Bg  are  derived  from  a 
single  potential  Z2-  Second,  the  coupling  term  between  field  and  plasma 
is  properly  defined  to  all  orders  in  the  radial  fluid  convection  velocity 
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g  «  Vf]uid^’  an<*  smooth1^  ac^'ts  marginally  stable  drift-speed  limited 

conduction  and  (quasi -static)  thermo-electric  contributions.  Properly,  as 

fast  time  scale  phenomena  which  track  the  instantaneous  state  of  the  system, 

these  effects  should  be  calculated  from  a  single  time  slice  of  information. 

Without  the  Hertz  vector  formulation,  an  electromagnetic  calculation  evolving 

E  .  separately  can  only  accomplish  this  to  the  extent  that  it  satisfies 

the  temporal  self-consistency  constraint  mentioned  above.  Third,  the  coupling 

term  Is  a  spatial  convolution  with  a  fixed  Green's  function.  The  spatial 

averaging  will  generally  serve  to  lower  the  noise  level  of  the  calculation 

by  damping  any  chatter  that  may  appear  on  the  discrete  source  array  extracted 

from  the  plasma  representation.  Fourth,  while  the  use  of  Z  does  require  the 

extraction  of  high  order  spatial  derivatives  from  discrete  data,  the  accuracy 

of  such  spatial  differencing  can  be  Improved  by  the  use  of  so-called  smooth 
9 

Interpolators.  The  estimation  of  local  E2  fields  from  the  parent  Hertz 
vector  by  mean:*  of  these  splines  is  routinely  accurate  to  [  10~61 ,  even  between 
the  input  data  points. 

The  Hertz-vector-based  electromagnetic  calculation  proceeds  from  a 

o 

spatially  discrete  set  of  Information  {Z,  3_Z,  J,  3TJ,  3TJ},  given  at  a  single 
time.  The  first  two  elements  of  this  set  are  defined  at  fixed  radial  points; 
the  last  three,  established  (sequentially)  on  a  convecting  (Lagrangian)  grid  of 
plasma  cell  centers.  The  local  values  of  E,  8  at  any  such  time  level  are  deter¬ 
mined  by  (Z,  3.Z,  J)  and  partially  determine  the  plasma  stress  and  heating  rates. 
From  the  values  of  3TJ,  3*J,  a  spatial  convolution  produces  the  Hertz  vector 
source  terms.  These  source  terms  are  interpolated  for  use  in  an  explicit 
wave  quadrature  formula  which  advances  the  Hertz  vector  and  its  partial  time 
derivative  in  order  to  establish  Z,  3TZ  at  any  subcycle  time  step  (used  in 
the  thermal  conduction  scheme)  or  at  the  next  major  time  level. 
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The  second  electrodynamic  option  is  that  of  E-field  diffusion,  which 
derives  from  the  Hertz  potential  wave  equation  in  the  limit  that  incoming  and 
outgoing  waves  are  in  detailed  balance.  This  limit  is  approximate  but 
reasonably  accurate.  It  neglects  the  net  displacement  current,  making  any 
detailed  outgoing  wave  calculation  impossible  because  no  wave  source  is 
extracted  from  the  plasma.  Reflected  power  is  manifested  solely  as  a 
diminished  voltage  on  the  plasma  load  and  the  energy  delivered  is  simply 
/tdt  t’plasna'W'  Here  0nly  the  plasma  current  waveform  is  available  as 
diagnostic  information,  insofar  as  the  voltage  at  the  convecting  plasma 
vacuum  Interface  is  not  usually  available  from  experiment.  The  major  simpli¬ 
fication  is  that  Z  need  not  be  explicitly  employed;  it  is  sufficient  to 
diffuse  the  E_  field  and  calculate  the  Bfl  field  in  a  (.self-consistent)  quasi- 

Z  U 

static  manner.  Diffusion  of  E„  rather  than  BQ  is  preferable  because  the 

z  y 

criteria  for  drift  speed  limitation  stem  most  directly  from  the  application  of 
an  "E-field-to-local-O"  mapping. 

Diffusion  based  electrodynamics  proceeds  concurrently  with  subcycled 

thermal  conduction.  A  circuit  relation  sets  the  time  varying  value  of  Ez  at 

the  outer  boundary  of  the  plasma,  and  the  time  slice  information  (Eth,  El  B, 

8,  a,  a}  defined  on  the  (convecting)  plasma  cell  centers  determines  the 

material  derivative  of  E„.  Here  E'  *  E  +  S„  x  Bfl  is  the  field  in  the 
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convecting  frame;  Eth  is  the  axial  thermoelectric  field.  The  new  field  Ez 
at  the  next  time  level  is  determined  implicitly  by  the  time  seauence  of 
material  derivatives.  The  plasma  mesh  evolves  in  response  to  the  stresses 
partially  determined  by  the  local  values  of  E,  B.  In  this  mode  of  operation 
some  iteration  is  required  in  obtaining  both  DE_/Dt  and  B,  in  contrast  to 
the  electromagnetic  calculation.  There  only  the  quasi-static  portion  of  B0 
requires  iterative  refinement  due  to  the  magnetic  field  dependence  of  the 
plasma  conductivity. 


In  either  electromagnetic  or  diffusive  mode  the  imposition  of  drift- 
speed  limited  conduction  is  straightforward  here  because  the  Ez(r,t)  profile 
is  established  a  priori,  before  J_  is  estimated.  The  rational  basis  for 
preference  between  these  two  can  be  derived  only  when  the  goals  of  I. A. (4) 
above  are  reached  and  when  the  comparative  computation  costs  of  well  tuned 
subroutines  are  established. 

The  hydrodynamic  problem  is  readily  formulated  (once  E,  B  are  con¬ 
sidered  given)  in  terms  of  well-known  techniques.  The  primary  challenge  is 
to  produce  a  quiet,  relatively  robust  algorithm  capable  of  concentrating  the 
computational  effort  on  the  faster  time  scales  of  the  problem,  minimizing  the 
number  of  expensive  derivative  evaluations  required,  and  retaining  spatial 
resolution  of  very  highly  compressed  plasma  loads. 

In  view  of  these  requirements  a  Lagrangian  formulation  of  the  two- 
fluid  transport  equations  for  number  density,  flow  velocity,  and  temperature 
has  been  implemented.  There  are  several  simplifications  to  the  complete  set 
of  transport  equations  which  are  appropriate  to  the  imploding  diode  plasma, 
but  the  hydrodynamic  model  is  essentially  that  of  Braginskii.  In  particu¬ 
lar  the  radial  flow  field  is  assumed  to  be  identical  for  ions  and  electrons, 


the  electron  density  is  assumed  to  be  in  quasineutral  equilibrium  with  the  ion 
density  (apart  from  small  ambipolar  charge  separations),  and  the  system  of 
moment  equations  is  closed  with  an  ideal  gas  equation  of  state  for  each 


fluid.  This  reduces  the  fluid  variables  to  (n^,  Vr,  Tj,  Tg,  s^}.  However 
Tg  is  not  a  very  good  thermodynamic  variable  because  it  is  tightly  coupled 
toy,  the  effective  ionization  state  and  to  £j,  the  specific  chemical 


potential.  In  order  to  follow  the  ionization  dynamics  efficiently,  let 


6fi  *  Te  TjT  and  *v0^ve  instead, 
radiative  equilibrium^  (CPE)  model  is 


The  approximation  of  the  collisional 
to  assume  a  very  rapid  establishment 
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of  ^and  from  Csf ven)  nj  and  T  .  Hence  the  usual  heating  sources  and 
sinks  for  T^,  when  ^is  fixed,  become  tin  effect)  heating  sources  and  sinks 
for  Sfi  in  the  context  of  CRE  ionization  dynamics.  From  a  sufficiently  accurate 
representation,  preprocessed  by  the  complete  CRE  model  and  smoothly  inter¬ 
polated,  one  may  simply  advance  0g  (.the  "grand  canonical  temperature")  and 
compute  Te,  £j  and  all  other  transport  coefficients  from  the  e-value 
obtained  at  any  timestep.  This  later  method  is  presently  implemented.  As  is 
common  practice  in  Lagrangian  calculations,  one  may  also  remove  the  ion 
density  from  the  set  of  evolved  quantities  byconserving  the  total  number  of 
Ions  in  any  cell,  unless  a  rezontng  is  done,  and  calculating  the  density  n^ 
from  the  time  varying  cell  volume.  Hence  the  irreducible  set  of  fluid  vari¬ 
ables  to  be  evolved  is  {Vr>  Tj,  eg},  with  the  r  =  defining  the 

variable  mesh  positions  to  be  used  in  computing  the  density  n^.  For  the 
choice  of  diffusive  electrodynamics,  this  set  is  expanded  to  {Ez>  Vr,  Tj  9e> 
and  the  separate  evolution  of  (Z,  3^  Z)  through  the  Hertz  vector  wave  equation 
is  discarded. 

The  complete  set  of  magnetoplasma,  ^-dependent  transport  coefficients 
Is  available  in  a  single  subroutine,  and  those  required  for  the  simplest 
relevant  problem  are  the  axial  electrical  conductivity,  the  radial  ion  and 
electron  thermal  conductivities  and  the  radial  thermoelectric  coefficient. 

The  expansion  of  the  model  to  include  radial  ambipolar  electric  fields, 
viscous  stresses  and  viscous  heating  only  requires  activation  of  the  axial 
friction  coefficient  and  five  viscosity  coefficients.  Moreover,  because  the 
relaxation  time  is  the  cornerstone  of  the  transport  calculation  an  extension 
into  more  strongly  coupled  plasma  domains  is  quite  straightforward  on  a 
cell-by-cell  basis.  The  use  of  a  modified  relaxation  time  based  on  the  local 
Coulomb  logarithm  is  a  useful  first  approximation  in  strongly  coupled  plasma 


transport^  and  can  be  relevant  here  in  some  simulations  of  wire  or  foil 
annular  implosions. 

The  hydrodynamic  calculation  across  a  major  time  step  proceeds 

from  the  time  slice  {NT,  V  .  Tt,  e  >.  Here  Nt  (the  time  invariant  total 

i  r  i  e  i 

number  of  cell  ions),  Tt  and  a  are  defined  at  the  cell  centers,  while  1/ 

I  e  r 

and  r(t)  are  defined  at  the  cell  interfaces.  From  this  basic  data  the  model 
calculates  Tg,  ej,  all  transport  coefficients,  all  heating  rates  and  all 
fluid  stress  terms  (after  E  and  B  have  been  established).  The  major  time 
step  is  selected  from  a  competition  among  the  magnetosonic  Courant  condition, 
fractional  cell  area  changes,  local  truncation  error  in  estimating  future 
accelerations  (OV^Dt),  and  the  most  rapidly  varying  electromagnetic  source 
terms.  Once  the  major  time  step  is  set, a  time  window  is  defined  over  which 
thermal  conduction,  compressional  heating,  radiative  heating  or  cooling,  and 
ion-electron  exchange  heating  are  subcycled  at  the  cell  centers  (under  an 
assumed  cell  boundary  evolution).  This  subcycle  is  Implemented  using  an 
(implicit)  variable  order  Adams  predictor-corrector  or  Gear  method.  At  the 
completion  of  the  subcycle  the  explicit  flow  field  advance  (and  implied  cell 
boundary  motion)  is  either  retained  or  iteratively  refined  (becoming  the 
starting  point  for  an  implicit  method).  If  desired  the  iteration  proceeds 
until  self-consistency  is  achieved.  The  self-consistency  criterion  is  solution 
of  the  first  order  (non-linear)  difference  equation  for  the  fluid  accelera¬ 
tion,  to  a  specified  tolerance  at  all  mesh  points.  For  the  diffusion  based 
electrodynamics  the  subcycle  includes  the  Ez  field  evolution  as  well. 

A  central  element  In  the  hydrodynamic  model  is  the  question  of 
radiation  transport.  The  algorithm  presently  implemented  is  a  compromise 

between  the  full  CRE  calculation  and  the  much  simplified  local  approxima- 
12 

tion  of  the  SPLAT  code.  First,  the  emission  function  preprocessed  by  the 
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CRE  model  is  represented  over  the  relevant  density  and  temperature  domain. 

One  such  emission  curve  is  required  for  each  radiation  category  one  wishes  to 
calculate.  A  simple  set  is:  {lines  hv  <  1  keV,  lines  hv  >  1  keV,  free- 
bound  continuum,  &  Bremsstrahlung}.  Next, for  each  radiation  category  a 
damping  coefficient  is  computed  for  each  cell.  These  damping  coefficients 
are  path  integrated  along  a  selected  ray  to  produce  a  probability  of  escape 
for  that  photon  category.  From  the  escape  probability  field  a  matrix  of 
coupling  coefficients  is  derived.  Once  the  coupling  matrix  is  given, 
selected  inner  products  with  the  vector  of  cell  emissions  for  a  particular 
category  produce  the  cell-to-cell  photon  exchange  and  the  net  radiative  loss 
to  the  plasma  from  any  zone  in  each  radiation  category. 

This  formulation  has  the  advantages  of  being  easily  expanded  into  a 
multiple  group  algorithm  and  easily  modified  to  include  other  elements.  At 
present  Ajl  and  Ar  are  available,  although  the  representation  of  0g  is  not 
yet  available  for  Ar.  The  disadvantages  lie  in  the  added  calculation  of  the 
coupling  matrix  and  in  the  need  to  establish  a  fresh  matrix  for  each  radiation 
category.  It  is  hoped  that  the  use  of  implicit  methods  for  the  overall  hydro- 
dynamic  advance  will  offset  the  cost  of  these  radiative  loss  calculations  by 
requiring  fewer  time  derivative  evaluations  to  achieve  useful  accuracy. 


II.  ACCURATE  METHOOS.  FOR  NUMERICAL  DIFFERENTIATION 
A.  Survey  of  the  Theory  of  Smooth  Interpol  ants 

All  Interpolation  methods  are  based  on  some  criterion  of  distance 
between  functions.  The  free  parameters  of  the  Interpolation  scheme  are  to  be 
chosen  In  such  a  way  as  to  mini mi ze  the  "distance"  between  the  data  function 
F(x)  and  the  Interpolant  z(x ) .  Smooth  interpolants  are  derived  by  Inventing 
a  distance  criterion  which  measures  the  oscillation  in  a  function  as  well  as 
Its  mean  value.  The  sequence  of  distance  measures  (or  norms) 

E  »  ^/b  dx  jF(x)  -  z(x)|2^  * 

Bn/  dx  |F(n)(x)  -  z(n)(x) l2^5  (II. I) 

Bn  /b  dx  wn(x)  |FCn)(x)  -  z(n>(x)|2)  15 

represents  a  progressively  more  general  set  of  functional  similarity  criteria. 
(Here  the  superscript  In  parentheses  indicates  the  ntb  derivative.)  Conven¬ 
tional  interpolation  schemes  stem  from  E  (.Hilbert's  norm)  while  smooth  inter¬ 
polants  arise  in  connection  with  T  (Talmi-Gilat  norm)  or  I.  to  be  discussed 

W 

here.  The  Bn,  wn(x)  are  rather  freely  chosen,  subject  only  to  the  constraint 
of  convergence. 

Each  of  these  norms  can  be  represented  also  as  an  inner  product 
operation  involving  smoothness  functionals,  i.e., 

! 

i 

1  I„  <9.  K)  I  E  B„  /b  dx  w„(x)  gCn)(x)  h(l,)(x),  (II. 2) 

n*  0  n  a  n 

i 

I  so  that,  e.g. ,  T2  5  1^  (F-z,  F-z).  It  is  these  smoothness  functionals  that 

provide  the  concrete  algorithms  to  be  used  for  interpolation.  In  particular. 
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choose  a  set  g^x)  of  complete  functions  (linearly  independent  and 
orthonormal  over  the  interval  [a,  b])*  Any  interpolant  z(x)  can  be 
expanded  in  terms  of  this  set,  and  the  expansion  coefficients  are  chosen  in 


such  a  way  as  to  minimize  a  smooth  norm  of  z(x)  (e.g.,  I  (z))  subject  to 

W 


constraints.  The  constraints  are  the  known  data  concerning  F(x), 
z(xj)  *  F(xj) ,  z  ^(x^)  *  F^(xk)  for  specified  sets  j  *  1,  2,  3...J, 
k  *  1,  2,  3...K,  n  *  0,  1,  2...N. 

This  defines  an  Euler-Lagrange  problem  which  can  be  transformed 
(using  the  orthonormal i ty  of  the  gk)  into  a  simple  linear  system  for  the 


Lagrange  multipliers  {X^}.  The  matrix  for  this  linear  system  is 


determined  by  evaluating 


£  9,(X)  g.(Y) 
R  (X,  Y)  2  E  i'  q  r 
a- o  lw  V 


(ii. 


and  Its  partial  derivatives  at  the  mesh  points,  i.e.,  x  *  {x^}, 

y  *  {x.}.  It  is  the  choice  for  B  and  w  (x),  therefore,  which  concretely 
j  n  n 

specifies  the  interpolation  spline  R(X,Y).  The  original  work  of  Talmi  and 

g 

Gilat  shows  several  sets  {B„}  which  are  useful  in  that  they  render  R(X,Y) 

n 

sunmable  in  closed  form.  A  more  versatile  result  is  obtained  when  one 
invokes  7  rather  than  only  7  as  they  did. 

W 

Instead  of  searching  for  a  set  {Bn>  giving  a  useful  R(X,Y),  reverse 
the  process  and  choose  a  spline  R(X,  Y)  *  R  (X-Y)  likely  to  be  similar  to 
the  smooth  functions  expected  to  support  the  data.  The  norm  ?w  is  easily 
shown  to  be  so  general  as  to  allow  the  {Bn>  to  be  extracted  from  the  spline 
R  (X-Y),  rather  than  the  original  (reversed)  procedure.  For  example,  select 
w  (x)  «  (1  -  x2)n,  [a,b]  *  [-1,  +1]  ,  and  g,(x)  *  PJl(x)  (a  Legendre 


Polynomial).  Then 


V9*’  gz}  * 

depends  only  on  those  Bm  forO  <  m  <  l,  Since  any  function  R(X,Y)  can 
be  expanded  on  [a,b]  in  terms  of  transformed  coordinates  U  =  X-Y,  V  =  X+Y 
one  may  always  adjust  the  {Bn}  such  that  the  g^X)  g^CY)  expansion  reproduces 
the  required  coefficients  for  the  expansion  of  R  in  terms  of  U  alone. 

Using  R(X,Y)  as  the  fundamental  element  one  may  cast  the  general 
spline  and  its  defining  linear  system  in  a  very  simple  form. 


Here  the  function  R(x,y)  is  assumed  available  in  closed  form,  the  F^ 
represent  input  data  giving  the  data  function  and  its  derivatives  (through 
order  o)  at  any  of  the  x^.  The  notation  convention  here  is  that  if  any 
Is  left  out,  so  also  is  its  corresponding  A^  . 

In  summary,  the  A^  represent  the  unknown  Lagrange  multipliers  re¬ 
quired  to  minimize  the  interpolant  with  respect  to  the  I  norm  implied  by  one's 

W 

choice  for  R(x,y),  subject  to  the  data  constraints  embodied  in  the  F^ml 

The  usual  application  of  this  involves  specifying  derivatives  only  at 
the  end  points,  but  this  is  only  a  matter  of  choice  and  one  may  readily  extend 
the  theory  to  include  integral  constraints  as  well  (m  <  0). 

A  compact  matrix  notation  for  this  algorithm  allows  one  to  state 
some  simple  rules  for  its  use  as  a  numerical  differentiation  scheme.  Let 


an; 


R(x. 


t 


4"1 

Rlx.xj)]  =  ^ 


» 


(generally  not  a  square  matrix), 


and  F^  5  F^,  xj^s  \  with  k*l,2,  3...  m  +  £,  . ..  J  +  N»(number  of 
derivative  conditions  given).  The  defining  linear  system  is  simply  R*X  *  F 
while  the  spline  at  any  point  x  becomes  z(x)  *  X  •  3^  R(x,y.).  Here 
the  operator  3^  Implies  a  derivative  only  for  those  elements  of  X 
corresponding  to  derivative  constraints  of  a  particular  order;  for  most 
X  elements  this  derivative  will  be  of  order  zero,  i.e.,  correspond  to  a 
function  value  constraint.  If  one  wishes  to  write  the  analog  of  a  finite 
difference  formula  estimating  the  derivatives  of  a  data  function  F  at  the 
mesh  points  used  in  R,  then  R-1  is  required. 


4P)  F  *  3iP)  1 


F 


F 


(II. 4) 


Similarly  the  method  provides  a  Jacobian  defining  the  sensitivity 
of  a  derivative  estimate  to  the  data  values  producing  it,  i.e. 


Derivative  estimates  off  the  original  data  points  are  given  by 
the  formula 
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(II. 6) 


In  comparison  with  more  common  interpolation  schemes,  the  "smooth" 
interpolators  have  several  advantages.  First,  the  algorithm  is  clearly 
vectorizable,  except  for  the  calculation  of  R-1.  In  a  piecewise  spline 
scheme  any  evaluation  point  x  (distinct  from  the  x.)  must  be  nested  in 

J 

its  appropriate  interval  before  the  interpolation  can  proceed  and  this 
represents  additional  (non-vectorizable)  calculation  overhead.  Using 
(II. 6)  to  evaluate  the  interpolant,  one  need  not  even  put  the  values  x  in 
any  particular  order.  A  further  advantage  over  piecewise  splines  is  a 
reduction  in  storage  space  for  the  spline  coefficients;  a  piecewise  scheme 
requires  the  storage  equivalent  of  M  «  (J  +  N),  where  M  is  the  order  of 
the  spline.  There  are  also  significant  increases  in  accuracy,  which  will 
be  discussed  below,  and  the  treatment  of  nonuniform  meshes  becomes  much 
less  formidable. 

B.  A  Family  of  Generalized  Splines  R(x,y) 

One  of  the  original  splines  noted  by  Talmi  and  Gilat  corresponds 
to  the  choice  Bn  3  Dn/n!  in  I,  and  leads  to  the  Gaussian  spline 

G(x,y)  «  — ^-j-  exp  -  (x-y)2/4D  .  (II. 7) 

("DP 

The  parameter  D  plays  the  role  of  a  correlation  width  in  coupling  the  data 
points  across  the  mesh,  and  thus  controls  the  condition  number  of  the 
resulting  matrix  Rg.  Qualitatively,  as  D  increases  from  some  small  value, 
the  Interpolant  Z(x)  evolves  from  a  "pickett  fence"  structure  toward  a 
smoother  structure,  then  further  increases  in  D  will  require  X  elements 
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w 


m 


m 


tty 


of  alternating  signs  in  order  to  follow  the  data  variations.  As  D  •*  » 
the  Rq  matrix  moves  toward  algorithmic  singularity  for  any  inversion  method 
using  finite  precision  arithmetic.  Usually,  for  positive  definite  data 
exhibiting  clear  trends,  an  optimal  0  will  be  found  (for  any  particular 
mesh)  which  allows  a  smooth  fit  with  positive  definite  values  X.  In  the 
case  of  non-uniform  meshes  it  is  useful  to  let  D  vary  across  the  mesh 
according  to  the  local  data  point  density.  The  spline  becomes 


G(x,  y„  0.)  ■  — - — r-  exp  -  (x-y.)2/4D, 

J  }  (trDjP  J  1 

so  that  the  effective  norm  now  involves  the  mean  value  F,  a  Bn  a,  IfVn!, 
and  some  weight  function  wn(x).  One  finds  in  general  that,  in  order  to 
construct  a  globally  smooth  interpolant  on  a  highly  distorted  mesh,  an 
adjustment  of  which  preserves  the  number  of  e-foldings  per  interval  is 
advisable.  Such  an  adjustment  acts  to  diminish  the  variations  in  the 
effective  bandwidth  of  RQ. 

A  few  numerical  examples  of  this  interpolation  method  with  the 
Gaussian  (or  G-spline)  are  sufficient  to  convince  the  user  that  functions 
F(x)  with  nearly  null  derivatives  (anywhere  on  the  domain)  are  rather  path¬ 
ological  cases.  This  leads  quite  naturally  to  the  exploration  of  other 
R(x,y)  in  attempting  to  accommodate  such  data.  A  useful  set  is  derived  by 
integrating  the  G-spline  repeatedly: 


E  (x,  yy  Dj)  *  1  +  erf  (u^) 


El(x>  *j*  V  *  2  °j  |uj  [  +  erf<uj>]  + 
E2(x,  yj,  Oj)  *  4Djj[l  +  erT(uj)j  +  ?]  + 


(II. 8) 
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where  u.  *  (x  -  y.)/2  0?  determines  the  relative  coupling  of  the  data 

J  J  J 

points.  The  E-splines  are  the  smooth-norm  generalizations  of  the  familiar 
8-spline  elements  (x  -  y)^  used  in  piecewise  spline  theory.  They  are  of 
particular  utility  in  accommodating  data  which  exhibit  a  clean  asymptotic 
trend  or  areas  of  null  derivatives.  For  example  in  Table  I  is  shown  a 
(fully  optimized  in  0)  test  of  the  E-spline  sampling  a  function 

/x  -  xA4 

F(x)  *  A  exp  -i — g - 1  +  C  exp  - 

at  a  progressively  larger  number  of  points.  The  quantity  tabulated  is  the 
average  fractional  error  in  for  m=0,  1,  2,  3,  A,  over  the  domain  [-1, 

+1 ].  The  interpolant  and  test  function  are  compared  on  a  mesh  much  finer  than 
that  used  to  sample  F  and  to  originate  the  calculation  of  X  values  for  Z(x). 


TABLE  I 


rel.  error  in  Z^, 

averaged  over  [-1,  +1]: 

Original  sample 
points  of  F 

m*0 

1  2 

3 

4 

21 

2.06E-3 

0.346  0.344 

0.58 

22.3 

31 

1.10E-5 

4.29E-4  1.16E-3 

8.26E-3 

3.96E-2 

41 

4.69E-7 

8.92E-5  1.92E-4 

6.79E-4 

2.30E-3 

51 

2.53E-7 

8.77E-6  6.09E-5 

7.16E-4 

7.77E-3 

61 

7.41E-7 

7.49E-5  2.59E-4 

3.15E-3 

8.25E-3 

As  is  evident  from  these  results  the  E-spline  can  track  accurately 
data  of  quite  different  functional  form.  This  example  is  also  notable  in  that 
a  G-spline  would  provide  a  much  poorer  result  -  the  constants  in  F(x)  had 
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been  adjusted  so  as  to  produce  a  very  flat  region  in  the  vicinity  of 
x  *  0.2.  The  degradation  in  this  interpolation  for  larger  numbers  of  sample 
points  is  a  manifestation  of  accumulated  roundoff  in  the  linear  system 
solution.  Moreover,  in  common  with  other  methods,  most  of  the  average  error 
is  coming  from  the  ends  of  the  domain.  The  accuracy  in  the  interior  is  often 
an  order  of  magnitude  better  than  the  mean. 

As  with  any  fresh  technique  there  are  some  points  less  well  under¬ 
stood  in  connection  with  these  new  splines.  In  particular,  one  is  a  robust 
algorithm  for  the  choice  of  D,.  At  present  the  complete  automation  of  these 
splines  is  not  quite  satisfactory.  Choosing  D.  based  on  the  sample  density 

J 

certainly  works  roughly  but  the  optimum  value  is  apparently  not  independent 
of  the  input  data  F,  and  more  systematic  study  will  be  required  on  this 
topic.  Also  the  E-splines  tend  to  overshoot  functions  which  rise  and  level 
off  sharply,  but  this  can  probably  be  corrected  by  adding  an  offset  in  the 
argument  u^,  i.e.  u^  «  (x  -  y  +  Oj)/2  Dj  . 

The  experience  thus  far  has  been  very  positive  with  respect  to 
"interactive"  fits  (i.e.  where  the  user  chooses  the  o.  and  D.  rather  than 

J  J 

automating  their  choice),  and  in  the  case  of  extracting  the  electric  field, 

E,  from  a  Hertz  potential  on  a  fixed  but  non-uniform  radial  mesh.  Specific 
examples  are  discussed  in  Section  D. 

C.  The  Interactive  Interpolation  Package  (VAX  Implementation) 

The  present  vehicle  for  preprocessing  CRE  emission  or  equation-of- 
state  data  is  a  pair  of  command  file  synonyms:  REPRESENT  and  EVALUATE.  The 
details  of  their  usage  are  documented  in  Appendix  I.  The  REPRESENT  algorithm 
processes  an  input  file  assumed  to  contain  a  mesh  y.  and  one  to  three 

J 

functions  F.  evaluated  at  the  mesh  points.  The  user  is  prompted  for  the  spline 

J 

choice  {G-splines,  E-splines,  or  E ^  splines}  and  the  width  parameter  D(or  D 
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if  the  given  mesh  is  found  to  be  non-uniform).  The  appropriate  R  matrix 

is  calculated  for  the  mesh  given  and  an  IMSl  linear  system  solver  CLEQT2F) 

is  used  to  calculate  a  A.  vector  for  each  function  F.  input.  The  output 

J  J 

consists  of  two  files  by  the  same  root  name  (user  supplied).  The  first  is 

suffixed  "H"  and  is-  a  header  file  providing  default  and  dimension  informa¬ 

tion  to  the  companion  process  EVALUATE.  The  second  is  an  interpolation  file 

containing  a  copy  of  the  original  mesh  y . ,  the  D.  required  and  the  resulting 

J 

Aj  for  each  of  the  functions  Fj  input. 

In  the  present  configuration  REPRESENT  allows  the  user  to  specify 

boundary  derivatives  of  first  and/or  second  order  at  either  limit  of  the  y 

domain.  The  boundary  condition  descriptors  are  stored  in  a  2  x  2  array,  the 

boundary  derivative  matrix  (80M),  and  embedded  in  the  header  file.  The  overall 

accuracy  of  the  linear  system  solution  is  controlled  by  a  specified  number 

of  digits  of  precision  in  the  input  file  (IDGT)  and  the  number  of  iterative 

refinements  allowed  to  the  IMSL  solver  (ITMAX) .  If  invoked  ,  the  iterative 

refinement  proceeds  until  the  solution  can  be  converged  to  IDGT  precision. 

Usually  iterative  refinement  is  not  needed  and,  for  the  G-spline,  often 

appears  to  be  unstable  as  well. 

The  companion  process  EVALUATE  produces  a  file  of  interpolant 

values  and  derivatives  through  some  specified  order  (.<_  4), for  any  of  the 

Aj  vectors  contained  in  an  output  file  from  REPRESENT.  The  evaluation  domain 

can  be  any  subset  of  the  original  mesh  Ly-i  <  y  <  y  )  or  can  even  line  outside 

'  J 

it.  EVALUATE  also  allows  the  user  to  change  the  spline  type  if  desired. 

For  example  A^  calculated  with  a  G-spline  can  be  evaluated  with  an  E-spline 
to  estimate  the  (definite)  integral  of  the  original  data  over  any  subinterval. 
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(i)  Estimating  an  Ez  field 


The  primary  problem  in  using  the  Hertz  potentia1  is  the  require¬ 
ment  that  higher  order  derivatives  be, extracted  from  the  discrete,  time- 
dependent  potential  Z ( r^ ,t).  As  a  test  of  the  smooth  interpolator's  ability 
to  resolve  these  fields,  a  model  Hertz  potential 


Fz<*>  -  r 


(1  +  e)x 


-/ 


d  u 


e"u-  1 


Jt  J* 


_  ri  ao 


where  x.  5  — p —  ,  was  chosen  as  a  data  function.  This  is  the  Hertz  potential 
1  b* 
for  the  E  field 


E(r,t)  ■  £Q(t)  |l  +  e  -  exp-(r2/b)J.  , 

exhibiting  a  sharp  rise  in  the  vicinity  of  r-b*5.  In  a  similar  manner  to 
the  test  described  in  Table  I  above,  the  function  Fz(x)  was  sampled  at  41 
points  on  a  uniform  x  mesh  (hence  a  mesh  in  r  which  is  denser  near  the 
origin).  The  interpolation  was  carried  out  using  the  E2  spline,  and  the 
interpolant  and  Its  derivatives  were  compared  to  those  of  the  original, 
analytically  given  test  function  on  a  mesh  much  finer  than  the  interpolation 
grid.  The  results  are  shown  in  Table  II  as  an  average  relative  error  and 
a  typical  interior  relative  error  for  the  interpolant  z(x)  and  its  first 
four  derivatives. 


TABLE  II 


2(nl)  m*0 

1 

2 

3 

4 

mean  error 

1.36E-5 

1.00E-5 

5.21E-5 

1.0E-3 

0.239 

interior  error 

2.47E-6 

3.25E-7 

2.73E-6 

7.56E-6 

1.03E-4 

Again  most  of  the  average  error  occurs  near  the  end  points.  The  test 
function  is  well  matched  even  between  the  sample  points  (where  one's 
Lagrangian  fluid  grid  will  almost  always  be  located). 

The  good  fit  to  the  potential  holds  as  well  in  estimates  of 
E(x)  and  3x  E,  even  though 

E  a  +  x  3*  zl  ,  and 

4a  XJ  l  x  x  ) 

o 

3*E  'IT7  (x23x2-x3xz-3v) 

o 

involve  linear  combinations  of  the  derivatives  of  Z.  Typical  interior 
errors  for  E  are  3E-6  and  for  3xE,  IE-6.  The  E2  spline  has  routinely 
extracted  these  E-fields  in  many  such  tests,  though  the  spline  function  by 
itself  only  roughly  resembles  the  test  function. 

(ii)  Representation  of  thermodynami c  data 

In  order  to  avoid  advancing  both  Tg  (electron  temperature)  and 

£j  (specific  chemical  potential)  concurrently  in  tin®,  the  use  of  a 

collisional  radiative  picture  will  hold  these  quantities  in  quasi -equilibrium. 

The  calculation  of  this  quasi -equilibrium  self-consistently  with  the  plasma 

evolution  is  a  very  expensive  process;  it  is  therefore  useful  to  fit  a 

paradigm  calculation  of  and  £  (the  mean  charge  state)  and  use  this  as  an 

approximate  representation  of  the  underlying  atomic  physics. 

The  fundamental  variable  used  is  a  "grand  canonical  temperature," 

2 

ee  2  Te  +  3  V*  It  can  be  most  easily  represented  by  defining  a  branching 

ratio  b  *  8  /T  and  examining  b(nT,9  ).  Owing  to  the  large  domain  of  9 
SC  IS  s 

values  relevant  to  the  model,  the  optimal  data  to  fit  are  In  b  and  in  9  . 


In  Figures  1  and  2  respectively  are  plotted  the  interpolants  to  in  b  vs. 

in  9g  (for  nj  *  10^ 9)  and  the  fits  to  £jCTg)  and  ^(Tg)  with  6g  shown  for 

19 

comparison  (again  at  n^  »  10  ) .  The  plots  are  very  true  to  the  original  dUa 

(5-6  significant  figures)  and  were  generated  by  splines  of  the  form 


in  b(9fi,nI)  »  £  xj  (nj) 

j 

?'-£n!?xf(ni) 

J 


G(£n  9  ,  In  9 . ,  D . ) 

"  J  J 


•  E(£n  Tg,  In  T.,  D) 


£j  *  -frij  2  X*  (nj) 


•  E (in  Tfi,  In  Tjt  D)] 


(II. 9a) 

(II. 9b) 

(II. 9c) 


with  a  sampling  of  31  data  points.  The  choice  of  G-splines  or  E-splines  is 
guided  by  the  general  form  of  the  data.  The  branching  ratio  b  is  a  cleanly 
peaked  function,  so  a  G-spline  is  the  natural  choice.  The  ionization  state 
and  chemical  potential  curves  exhibit  fla"  oortions  (in  between  the  opening 
of  atomic  shells)  so  the  E-splines  are  indicated. 


III.  THE  HERTZ  DRIVEN  IMPLODING  PLASMA  RADIATOR  (PROGRAM  ZDIPR) 


A.  Overview 

A  computer  program  is  being  constructed  to  solve  the  one-dimensional 
(in  radius)  time-dependent  interaction  of  a  plasma  annulus  (with  temperature- 
dependent  degree  of  ionization)  and  an  electromagnetic  driving  field  in  a  cy¬ 
lindrical  cavity.  It  was  shown  previously  that  the  number  of  equations  describ¬ 
ing  the  fields  can  be  reduced  and  the  computational  stability  probably  improved 
by  the  use  of  a  generalized  Hertz  vector,  Z,  from  which  E  and  8  field  vectors 
can  be  derived  by  differentiations.  In  the  limit  of  high  plasma  density  the 
wavelike  term  in  the  field  equations  can  be  neglected;  this  corresponds  to 
a  field-diffusion  approximation.  In  simple  cases  —  constant  scalar  con¬ 
ductivity  for  example— this  gives  the  usual  magnetic  diffusion  equation. 

The  program  must  couple  field  or  Hertz  potential  quantities, 
evaluated  most  conveniently  on  a  stationary  (Eulerian)  grid  of  points,  with 
fluid  equations  describing  the  mass-motion  (V)  electron  and  ion  temperatures 
(Tg,  Ij),  ionization  (y) ,  and  effective  conductivity  (I)  of  the  plasma,  all 
functions  of  r  and  t  evaluated  most  conveniently  on  a  comoving  CLagrangian) 
grid.  The  differential  equations  are  Integrated  forward  in  time  by  a 
varlable-timestep  method  (GEAR),,  but  because  the  plasma  thermodynamic  state 
evolves  much  more  rapidly  than  the  fields  do,  it  is  updated  separately,  on  a 
fast  timescale,  with  the  subcycle  timesteps  shorter  than  the  major  timestep  T 
required  for  Integration  of  the  field  and  momentum  equations.  Doing  these 
updates  on  different  timescales  is  expected  to  save  a  significant  amount  of 
running  time  and  expense  in  the  computation,  and  may  make  the  difference 
between  a  practical  cost  and  an  impractically  large  one. 

In  the  code  description  that  follows, the  subroutines  that  update 
the  plasma  state  on  the  Lagrangian  "r"  mesh  are  called  MESHSTRESS 
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(forces) ,TDOT  (heating),  CREOS  (ionization  state)  and  FLUIDOTS  (material 
derivatives),  all  of  these  being  used  by  the  subroutine  GEARBOX  which 
controls  the  updating  as  part  of  the  overall  fluid-advance  routine  HYDROPUSH. 

The  subroutine  that  updates  the  field  or  Hertz  potential  variables 
on  the  Eulerian  "R"  mesh  is  called  ZPUSH.  When  the  Hertz  potential  represent 
atlon  is  used,  the  subroutine  HERTZDOTS  updates  the  source  terms  in  the 
potential  equation  at  each  major  timestep.  When  the  electric  field  diffusion 
representation  is  used,  the  corresponding  subroutine  EDOT  updates  the  field 
source  3tJ  and  alternate  plasma  state  movers  (.GEARBOX  .►  TETGEAR;  FLUIDOTS  -*■ 
TETDOTS)  are  used  because  different  electromagnetic  information  is  required. 

This  code,  called  ZDIPR,  is  now  mostly  written  and  some  of  its 
portions  have  been  separately  checked.  Further  testing  and  integration 
of  the  various  subroutines  will  be  carried  on  in  the  coming  weeks. 


B.  Theoretical  Summary 

The  electromagnetic  calculations  required  are  either  solutions  of 
the  Hertz  wave  equation  or  of  the  electrodiffusive  equation  derived  from 
it.  If  a  complete  calculation  is  desired,  all  fields  and  currents  are  derived 
from  solutions  to 


3?z  -  x'1  3U  (“  3U  Z)*  (y  3t  J  , 


(III.1) 


where  t  3  ct/R  .  u  *  R/Rrt,  Z  3  ra/ c,  J  3  IE,  E  3  E,/QaR"2,  B  3  Ba/Qr~2, 

O  w  0  ZOO  uQ 

£  *  -u  1  3u(u3u  Z)and  Z  3  Z2(r,t)/QQ  define  the  dimensionless  variables  and 
fields;  If  one  gives  up  some  information  concerning  the  details  of  the 
diode  fields  and  makes  the  assumption  that  the  incoming  and  outgoing  wave 
components  are  in. detailed  balance,  then  the  wave  equation  above  can  be 
transformed  to  a  diffusion  eouation 


h  E  ■  {u'1  3u  u3u  E>  -  E  3t  **  Z  -  3t  Eth  -  8  3t6  •  <In-2> 


where  Eth  Is  the  dimensionless  thermoelectric  field,  8  3  V^-|u1-d/c,  E'=  E  +  SB, 
*  D 

E  3  E  +  88  +  E^  and  5^  +  83u.  (This  diffusive  limit  is  derived  in 
Appendix  II.) 

The  fluid  response  to  the  electromagnetic  stresses  and  heating  is 
embodied  in  the  relations 


jjrV-Z  Kin  n.  - 

m  r  1  m  men  j  mn  j 


(III. 3) 
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3jnx  v°z  z‘ 


-  u 

mI  \  Te  / 


{III. 4a) 
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(III. 4b) 


In  these  expressions  eg  =  Tg  +  >  T  *  Tj  +  ^Te,  m  *  mg; 

X.  t  is  the  thermal  conductivity,  t  the  plasma  relaxation  time,  E„  the 
ambipolar  radial  field  (with  p  its  induced  charge  density),  and  Qra-d  is  the 
net  (local)  radiative  heating  or  cooling.  The  dimensional  version  of  fields 
are  subscripted  with  a  vector  component;  dimensionless  fields  are  not.  The 
radial  electric  field  is  established  as  a  solution  to  an  integral  equation 
derived  from  the  radial  component  of  Ohm's  law,  (cf.  reference  6).  The 
drift-speed-limited  current  condition  is  supplemented  by  a  (nonlinear) 
change  in  E  where  the  local  E  field  requires  it.  The  overall  architecture 
of  the  model  is  indicated  in  the  flow  chart  of  Fiqure  3,  and  some  detailed 
documentation  of  the  physical  content  has  been  mentioned  above  (Chapter  I). 

The  following  sections  describe  in  some  detail  all  the  specific 
algorithms  which  are  required  to  implement  CHI .1-4).  For  the  sake  of 
simplicity,  a  common  notation  for  all  of  them  is  given  here. 

The  (nonuniform)  time  levels  for  any  variable  are  indexed  with  a 
leading  superscript;  the  lagrangian  fluid  mesh  is  denoted  by  r  and  it. 
material  derivatives,  by  a  superscripted  dot  (or  dots).  Spatial  indexing  is 
denoted  by  a  trailing  subscript  and  various  cell-to-cell  averaging  operations 
are  denoted  by  an  overbar  or  by  angle  brackets. 
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Since  the  code  performs  a  subcycle  calculation  wherein  the  time 


step  is  set  by  the  IMSL  Gear  integrator,  these  intermediate  time  levels  are 
denoted  by  a  leading  superscript  asterisk.  The  subcycle  interval  is  denoted 
by  T,  and  its  end  points  by  Tlo  and  Thi.  Other  spatial  grids  are  employed 
in  addition  to  r^.  The  Eulerian  electromagnetic  grid  is  denoted  by  R^,  and 
the  (Lagrangian)  area  coordinate  a=r  or  (Eulerian)  coordinate  y=i/ua(R/Ro) 2 
will  be  used  also.  The  evaluation  of  any  quantity  at  the  Lagrangian  cell 
center  will  be  denoted  by  a  leading  subscript  ,rc".  For  example, *fj  denotes  a 
comoving  acceleration  on  the  boundary  of  the  cell  evaluated  at  some 
arbitrary  point  in  the  subcycle:  Tlo  <  x*  <  Tlo  +  T=Thi ;  and  13tZj(  denotes 
the  first  (partial)  t  -  derivative  of  the  Hertz  potential  at  the  i  major 
time  level  and  the  kth  grid  point  of  the  R  mesh.  On  the  other  hand, 

denotes  the  first  material  derivative  of  Ez  at  the  jth  cell  center  at  some 
time  t*  in  the  subcycle.  In  general  the  fixed  mesh  (R  or  u  or  y)  will  be 
indexed  spatially  by  k;  the  comoving  mesh  (ror  a),  by  j.  Exceptions  will 
be  noted. 


C.  Electromagnetic  Algorithms  and  Subroutines 

Either  the  electromagnetic  mode  or  the  electrodiffusive  mode 
requires  the  establishment  of  a  local  Ez  field  on  the  plasma  mesh  continuously 
through  the  subcycle.  The  electrodiffusive  mode  is  discussed  in  Section 
III.G(iv)  below;  only  the  electromagnetic  calculation  requires  distinct 
subroutines  and  is  discussed  here. 

The  intermediate  time  t*  chosen  by  the  subcycle  integrator  is  the 
fundamental  parameter  controlling  the  Hertz  potential  wave  equation  integrator 
ZPUSH.  The  input  data  to  ZPUSH  are,  in  addition  to  the  previous  Hertz 
potential  (with  its  partial  r  derivative)  at  the  last  major  time  level  and 
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the  source  terms  arising  on  the  Lagrangian  plasma  mesh,  i.e.  {1Zk>  ''a  Z^, 

13TCJj.  The  output  of  ZPUSH  is  simply  * 1 k,  *3.^,  from  which  one 

may  derive  the  required  local  fields  in  the  fluid. 

The  value  of  t*  determines  a  radius  of  integration  for  the  wave 
quadrature  formulae/  which  generate  the  time-advanced  potential  and  its 
time  derivative.  The  source  terms  and  initial  data  are  functions  of  radius 
only,  but  the  wave  equation  solution  requires  them  to  be  mapped  onto  a  regular 

array  of  quadrature  points.  The  quadrature  points  are  arranged  within  circles 

of  raoius  t*  centered  at  each  electromagnetic  mesh  point  R^.  The  source  terms 
from  the  plasma  mesh  V.  are  smoothly  interpolated  to  the  whole  space  R  using 

J 

a  G-spline  and,  when  combined  with  the  existing  interpolation  coefficients  for 
the  initial  condition  data,  one  may  define  six  wave  quadrature  functions, 

QF1 .... QF6 .  These  are  then  input  to  WAVEQUAD,  which  evaluates  them  at  the 
needed  quadrature  points  and  forms  the  inner  products  that  estimate  the  Poisson 
integral  solution  for  the  Z  wave  equation. 

Once  *Zk>  *3rZk  are  output  by  ZPUSH,  the  companion  routines  EFIELD 
and  8FIELQ  use  these  values  and  *J .  to  produce  *E.  and  *B.  by  means  of  E- 

J  J  J  • 

spline  interpolation.  Since  *J .  is  only  known  from  *E.  and  *1 .  the  previous 

J  J  J 

value  for  S.  is  used  to  estimate  *£ .  first.  The  quasi-static  B.  component 

J  J  J 

is  then  iterated  to  convergence,  as  discussed  in  Section  III.G(ii),  and 
a  self-consistent  *Ej  and  *8^  a.e  established. 

Ouring  a  subcycle  the  use  of  ZPUSH,  EFIELD  and  BFIELD,  with  , 

iC  j 

i  0 

®TCJj  fixed, provides  all  the  field  information  needed  by  the  hydro- 
dynamic  portions  of  the  program.  These  sources  are  updated  however  at  each 
major  time  level  i,  i  +  1,  i  +2...  using  the  subroutine  HERTZDOTS,  the 
fluid-to-field  interface  code.  HERTZDOTS  references  the  material  derivatives 
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1  •  f  ••  i  i  • 

J.  and  J.  and  the  plasma  mesh  r. ,  r.  in  order  to  compute  these  sources. 

J  CJ  C  J  C  J 

The  higher  order  material  derivatives  in  J  are  calculated  using  complete 
information  on  the  time  derivatives  of  E,  Z  and  the  drift  speed  limitations 
to  J. 

D-  Transport  Coefficients 

The  central  elements  of  the  field-to-fluid  interface  are  the 
transport  coefficients,  calculated  by  the  subroutine  BRAGINSKII.  This  code 
produces  values  for  as  many-  as  20  transport  coefficients,  as  functions  of  the 
ion  density  Cn j ) »  electron  temperature  (T  ) ,  ion  temperature  (T^ ) ,  an  effective 

c  field  Bq.  Input  units  are  CGS, 

except  for  temperatures  in  [keV]. 

The  coefficients  are  evaluated  at  the  cell  centers  and  a  local 
coordinate  system:  z  *  fQ/iB0|,  x  -  any  orthogonal  direction,  y  *  z  x  x  is 
employed.  The  z  is  called  the  parallel  direction;  x,  the  perpendicular;  and 
y,  the  cross  product.  Output  [CGS]  transport  coefficients  are  written  in  three 
common  blocks  (PARATRAN,  PERPTRAN,  CRPXTRAN)  which  connote  these  directions, 
and  in  two  additional  conmon  blocks  (VISCTRAN,  EXCHANGE)  which  contain  viscosity 
and  electron/ion  heat  exchange  rates.  The  logical  matrix  of  coefficient 
selectors  CONTROL  (6,  3)  in  the  common  block  TRANSPORT-CHOICE  determines 
which  coefficients  are  evaluated,  as  shown  below,  in  Braginskii's  notation. 


ionization  state  (j),  and  the  local  megneti 
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TABLE  HI 


Coefficient 

0,1) 

0,  2) 

0,  3) 

of  friction 

Ct„  (J,  B) 

CLjJ,  8) 

“a<?-  b> 

Thermoelectric 

(2,  1) 

(2,  2) 

C2,  3) 

coefficient 

sf  tj,  B) 

i'f  ij,  B) 

ba  9-  81 

Electron  thermal 

(3,  1) 

(3,  2) 

(3,  3) 

conduction 

X®  (J.  B) 

X®  (J.  B) 

v  9- 81 

Ion  thermal 

(4,  1) 

(4,  2) 

(4,  3) 

conducti on 

xl  (J.  B) 

x[  (J.  B) 

co 

-  < 

X 

Electrical 

(5,  1) 

(.5,  2) 

(5,  3) 

conducti vi  ty 

a,,  (j,  B) 

aA^B) 

Viscosity 

(6,  1) 

(6,  2) 

(6,  3) 

(ion  &  electron) 

V4  8> 

v«  81 

V4  9- 81 

If  CONTROL  (I,  J)  has  the  value  "true",  this  triggers  the  evaluation  of  the 

indicated  coefficient  at  all  cell  centers  of  the  input  mesh.  Any  true  value 

in  CONTROL  (6,  1,  2  or  3)  causes  evaluation  of  all  viscosity  coefficients, 

and  the  exchange  heating  is  always  evaluated  (independent  of  CONTROL). 

Since  Braginskii's  theory13  provides  all  these  coefficients  as 

functions  of  wjTp  weTg  and  the  evaluation  scheme  also  provides  arrays 

in  common  blocks  which  show  u>  t  ,  u )Txr,  t  .  rT  and  A,  the  Coulomb  logarithm, 

e  e  i  i  e  i 

at  each  cell  center.  The ^ dependence  tabulated  by  Braginskii  has  been 
parameterized  by  power  laws  in  2.’*.  These  fits  are  implemented  in  the 


35 


companion  routines  BRALPHA,  BRBETA,  BRGAMMA,  and  BRDELTA,  evaluating 

<4  <?).  ao  (?>• <?»•  (;)■  <;>• s;  <;»•  <;>•  <?). 

Yq  (^)h  and{5Q{^),  (^)>  respectively.  The  constant  values  a!f ,  6^', 

Y^'  are  also  returned  in  the  various  output  arrays.  All  of  these  ^dependent 

functions  can  be  represented  to  reasonably  accuracy  (^2%)  by  expressions  of 

the  form  cQ  +  c-j  (^)~*  on  [l,  +  «]  in  For  the  purposes  of  transport 

coefficient  evaluation  ^  is  restricted  to  be  greater  than  1  and  so  also  is 

the  Coulomb  logarithm.  However  the  systematic  modification  of  t  and  tT 

e  l 

allows  the  extension  of  the  theory  and  the  code  into  domains  of  stronger 

3 

coupling  (in  the  plasma  parameter  nXQ).  Once  the  appropriate  relaxation 
times  are  available  the  restriction  on  the  Coulomb  logarithm  can  be  relaxed. 

The  modification  of  tgdueto  turbulence  in  the  low  density,  drift 
speed  limited  regions  is  also  a  natural  temptation.  It  represents  the 
simplest  way  of  obtaining  a  consistent  and  systematic  transport  package  for 
turbulent  loads.  At  present  it  has  not  been  done  because  U)  the  model  of 
turbulent  relaxation  by  a  Fokker-Planck  kinetic  equation  is  questionable  on 
many  grounds,  (. U )  the  magnetic  field  effects  in  turbulent  relaxation  are 
not  likely  to  be  accurately  represented  by  Braginskii's  functions,  and 
(jLLL)  the  proper  treatment  of  turbulent  transport  will  depend  on  isolating 
the  detailed  properties  of  those  microinstabilities  peculiar  to  the  cylindri 
cal,  highly-inhomogeneous  plasma.  Those  properties,  which  can  be  calculated 
only  when  the  background  quasi-equilibria  are  available,  depend  on  the 
results  of  the  present  model  and  cannot  be  specified  a  priori.  The  marginal 
stability  criterion  J  <_  eng  cg  circumvents  these  difficulties  and  is  thus 
used  in  place  of  xa  modifications. 


E.  fluid  Oensitv  and  Flow  Field  Time  Derivatives 

The  data  base  for  all  the  hydrodynamic  calculations  is  the  correnon 

block  FLUIO-STATE,  containing  {N.,  V  . ,  V. ,  V.}.  An  update  of 

jci*jce,j  j  j 

the  FLUID-STATE  is  the  central  result  of  a  major  time  step.  The  relation¬ 
ships  among  the  basic  fluid  variables  on  the  lagrangian  mesh  are  illustrated 
in  Figure  4.  The  simple  two  and  three  point  area-weighted  differencing 
schemes  for  such  a  mesh  are  discussed  in  Appendix  III.  The  cell  ions 
{N - >  are  a  conserved  vector  of  ions/cm  resident  in  the  (compressible)  cell 

J 

^r^,  Vj+1],  assuring  strict  particle  conservation  and  a  solution  of  the 
equation  of  continuity  limited  in  accuracy  only  by  the  evolved  values  of 
{  r^}.  The  elements  of  {N^>  are  assigned  soatial  locations  given  by  the  cell 


{  r^}.  The  elements  of  {N^>  are  assigned  soatial  locations  given  b 
center  position  (defined  by  the  equal  area  point)  *r^  *  ’ 
They  change  only  if  a  regridding  Is  called  for. 


rj+l^ 


Boundary  Variables: 


Cantered  Variables:  QN. 

cT.j  cJi  A« 

p  □  Transport 


.  e®«j  cZi 


rs  :  ri+i 


i  1  1  1  1  1  1  1  1  1  1  ! 

cp  O  O  O  O  O  O  O  O  O  O  O  O  O  O  C)  C)  ffi  ia 

!  ! 

I  I 

S _ h  ai!ai+i  !'  'I  I" :  J  k  !  ! 


h<  •  h> 

Mesh 

Spacing 

Notation 


Figure  4.  The  Lagrangian  Fluid  Mesh 
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When  needed,  the  stresses  at  the  cell  interfaces  are  calculated  by 
the  subroutine  MESHTRESS.  This  code  has  two  modes  of  operation.  In  the 
usual  mode,  during  a  subcycle, the  acceleration  of  the  interface  *r .  is 

J 

computed.  But  at  a  major  time  step,  the  initiation  of  the  subcycle  and  the 
calculation  of  the  Z  field  sources  requires  an  explicit  calculation  of  the 
jerk  Vj  as  well.  The  fluid  stress  has  three  major  components  as  noted  in  Ea 
(III. 3)  above.  Each  of  these  stress  components:  pressure  gradient,  magneto- 
strictive,  and  ambipolar  electrostatic,  has  a  particular  realization  in 
terms  of  area  weighted  finite  differences  (or  smooth  interpolants  if  more 
precision  is  needed)  that  is  best  suited  to  the  mathematical  and  programming 
requirements  of  the  problem. 


First,  the  pressure  gradient  stresses  are  most  naturally  derived 


from  the  mesh  structure  shown  in  Figure  4. 
mean  fluid  density  at  the  cell  center  ^n^ 


Using  the  {N.}  and  {^a.}  the 

J  J 

■  v(  Vi  •  ad  )tt  is  the  quantity 


determining  the  local  inertia  of  the  fluid.  The  pressure  gradient  stress 

2 

is  thus  easily  shown  to  be  estimated  (to  0(Ar  ))  by 


j  'fluid 


2  r.  J- 
J  m 


caj  ‘  caj-l 


(III. 5) 


-  (2  r^/m) 


where  the  area-weighted  interface  temperature  T.  is  defined  by  (cf. 

J 

Appendix  III  also)  : 


+  h_:_c*izL)  t.  \ 

Rcaj  '  caj-l/  C  \ca3  '  caj-l/  C  J  ] 


(III. 6) 


for  j  on  the  index  domain  [2,  J]  with  J  the  maximum  number  of  cells.  If 
one  Insists, as  in  Section  III. G, on  a  null  heat  flow  across  the  outermost 
boundaries  r^  and  rj+j,  then  the  natural  boundary  condition  on  this  pressure 
gradient  stress  is  to  neglect  the  temperature  difference  term  and  simply 
extrapolate  the  logarithmic  density  gradient  from  boundary  2  and  boundary  J. 
The  boundary  pressure  gradient  stresses  are  then  estimated  by 


1  'fluid 


"12 'U  "III  j 

lca2  '  cal  i  c  1 


(III. 7a) 


tcaJ  “  caJ-l  ) 


?J+1*  fluid  *  “2  rj+l 


cT>  • 


(III. 7b) 


Next  one  must  estimate  the  magnetostrictive  stresses  using  J. 

~  3 

and  c8j  and  suitable  area  weightings.  The  same  averaging  operation  shown  in 
(II 1.6)  provides  the  formulation 


‘j'jB  •  -(^) 


(III. 8a) 


with  the  boundary  conditions 


a  I  ,  (sjzdl  CBe,l\ 

"i,jb  to.  1  J 

r  I  *  . (--id - c  ° 

rJ+l ■  JB  \c  rcnI(J  J 


(III. 8b) 


(III  .8C) 
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The  ambipolar  stresses  are  somewhat  more  delicate  in  that  their 
sources  are  drawn  from  both  zone  boundaries  and  zone  centers.  The  strict 
application  of  the  radial  Ohm's  law  results  in  an  integral  equation  for 
the  relative  charge  separation  at  any  spatial  location.  To  a  good  approxi¬ 
mation  however  one  may  neglect  the  recursive  nature  of  that  relation  and 
simply  evaluate  the  radial  component  of  E  at  the  cell  boundaries  in  a  quasi¬ 
static  sense  at  any  time  level.  In  terms  of  the  axial  electron  drift  speed 
u2  and  the  (dimensionless)  Braginskii  transport  coefficients  aA  this 
electric  field  is 


iHr1 4  {[•,(> 

■  sf  TI  3r  &  "I  +  Te  V*"  ?nl]} 


’  *7 Tj) 


(III. 9) 


and  can  be  defined  on  the  cell  boundaries  easily,  differencing  for  the 


gradients  in  n^,  nj  and  T  -  ~  Tj  to  the  boundaries  directly  and  averaging  the 
strictly  central  terms  to  the  intermediate  boundary.  The  boundary  con¬ 
ditions  on  E,  are  similar  to  those  on  the  maonetostrictive  stresses  with  the 
jr 

addition  that,  in  compatibility  with  the  thermal  conduction  algorithm 
discussed  below,  the  temperature  gradient  source  terms  are  assumed  to 


vanish. 


The  ambipolar  space  charge  is  hence  p*  jr  7  •  E  with  E  given 


on  the  cell  boundaries.  But  p  is  most  naturally  computed  at  the  cell 
centers  by  means  of 


(III. 10) 
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fluid 


+  V 


but,  at  the  major  time  levels  when  the  jerk  is  needed  the  subroutine  goes 
further  and  computes 


In  order  to  compute  the  heating  rates  at  the  cell  centers  from  the 
FLUID-STATE,  several  distinct  steps  must  be  taken.  These  distinct  calcula¬ 
tions  are,  however,  collected  in  the  single  subroutine  TDOT.  First  one  must 
recover  the  electron  temperature,  mean  ionization  state  and  chemical 
potential.  This  is  done  by  the  subroutine  CREOS  (Collisional  Radiative 
Equation  Of  State)  to  which  is  input  the  ion  density  (*n,  .)  and  the 

C  1  *j 

electron  “grand  canonical  temperature"  (*6  .).  Using  the  smooth  inter- 

C  6  9  J 

polants  discussed  in  Section  (II.D.ii),  CREOS  first  obtains  the  logarithm 
of  the  branching  ratio  and  calculates  using  (II. 9a) 

In  (JTe>j)  -  t*  (*0efJ)  -  In  Cjbj),  (III. 12) 

from  which  *T  .  is  obtained  through  exponentiation.  Both  *2.  and 

C  6fJ  C/j 

*CT  _•  are  then  obtained  using  (II.9b,c)  with  the  logarithmic  argument 

implied  by  the  relation  (III. 12).  In  addition,  CREOS  calculates  the 

derivatives  of  the  branching  ratio  with  respect  to  0g  and  Hj  so  that  one 

may  calculate  material  derivatives  of  9  and  T  using 

e  e 


The  material  derivatives  of  Tg  are  needed  for  the  evaluation  of  the  jerk 
as  well  as  for  diagnostic  purposes.  If  the  jerk  evaluation  is  being  done, 


2  2 

then  CREOS  continues  with  a  calculation  of  3  In  b/3(£n  9  )  and 
32  In  b /3(£n  nj)2  as  well. 

The  next  operation  is  the  establishment  of  the  transport  coeffic¬ 
ients  on  the  mesh  using  the  data  base  {*Bfl  *nT  *Tt  . ,  *T  . , 

*^j}  supplied  in  part  by  CREOS.  This  procedure  is  complicated  by  the  need 
to  iterate  on  *BQj  (through  *Jj)  because  the  conductivity  depends  on 
the  magnetic  field.  The  details  of  this  iterative  refinement  are  discussed 
in  Section  III.G(ii)  in  the  appropriate  context  of  the  subcycle  algorithm. 
In  essence,  any  call  to  BRAGINSKII,  given  the  data  base  defined  above,  will 
produce  the  needed  transport  coefficients  to  go  further. 

The  establishment  of  the  transport  coefficients  and  temperatures 
is  necessary  to  calculate  the  thermal  diffusion  rates.  This  calculation  is 


presently  done  in  conjunction  with  the  compression  heating  ( *Q 


(pdV)j 


-*(V  •  r)j  JTj)  by  the  subroutines  VTGRADS  and  VTDIFDER.  The  module  VT GRACE 
(Velocity  and  Thermal  GRADient  Sources)  produces  arrays  containing  the  area 
meshes ?<ly  and  the  discrete  velocity,  temperature  and  conductivity 

data  {£Tej»  J^ej*  £Xij»  cTIj*  *n  a  su1table  form  for  use  wltb  either 
the  smooth  interpolation  package  or  the  area-weighted  differencing  module 
VTDIFDER  (Velocity  and  Thermal  OIFferenced  DERi vati ves ) .  The  boundary 
conditions  enforced  by  VTDIFDER  are  those  of  null  heat  flow  across  the  outer 
boundaries  *r^  and  *rJ+1.  This  is  done  by  means  of  ghost  points  playing  the 
role  of  cell  centers  *r  ,  *r..1  (cf.  Fig.  4)  and  the  use  of  simple  variable- 

CO  C 

mesh  three-point  differencing  formulae.  This  yields  the  output  array 

TOOTGF:  T.lj.  JO, '.T,),.  JO*), .  JO*,),.  JO.TjJj.  JO*  Tjlj. 

*(3.  a)j}-  These  results  are  used  in  calculating  the  thermal  diffusion  rates 
c  a  j 


*(Qe>j  "  *(4  £  3a  Te  +  4a^3a  9a  Te  +  >1  V>j 


(III. 14a) 
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(111.14b) 


'(Qr)j  -  *(4  3a  Te  *  4a  <8ax]»a  VxJ^TjJJj 


The  compression  rate  (7  • V  )  is  of  course  the  last  item  listed  above  in 
TDOTGF. 

With  the  gradients  available,  one  may  now  evaluate  the  thermoelectric 
fields  and  calculate  a  local  drift  speed  for  the  conduction  current  Jz 
using  the  free  electron  density,  the  complete  E  field  E2’EZ  +  Etl,'  and 
the  classical  conductivity  a  .  At  each  cell  center  this  drift  speed  is 
compared  to  the  local  sound  speed;  if  greater,  it  is  lowered  to  the  marginal 

A, 

stability  limit  uz^cs*  With  Ez  fixed,  the  drift-speed  limit  thus  imposes 
a  decrease  in  the  effective  conductivity  of  any  cell.  This  new  conductivity 

A 

is  set  such  that  E,  •  *  en  c  ,  and  it  replaces  the  classical  conduct - 

z  turb  e  s 

ivity  orginally  computed.  The  ohmic  heating  is  also  radically  altered  when 
the  drift  speed  limit  Is  Imposed.  It  is  calculated  here  by  means  of  the 
product  Jz  •  Ez',  which  is  explicitly  indpendent  of  the  conductivity  but 
reflects  It  in  any  limitations  on  Jz-  The  usual  heating  rates  are  thus 
completely  specified  by  the  calculations  sumnarized  thus  far;  only  the 
radiative  coupling  remains. 

The  radiative  energy  transport  calculation  must  be  done  for  each 
radiative  category  one  wishes  to  treat.  At  present  a  vector  of  emission 
rates  is  calculated  by  the  function  CREMIT  for  each  of  four  classifica¬ 
tions:  line  radiation  hv  <1  keV,  line  radiation  hv  >1  keV,  free-bound 
continuum,  free-free  continuum.  This  function  contains  preprocessed  fits  to 
the  CRE  model  emission  strengths  in  lines  and  free-bound  continuum,  and  a 
simple  free-free  emission  model.  A  companion  function  DAMPIT  calculates  a 
vector  of  damping  rates  *d^  using  attenuation  estimates  appropriate  to  each 
radiation  category.  Once  the  data  base  {*e^,  *d^ ,  *r \ ,  £10.=  *CL+1-*G.  } 
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is  computed  for  a  photon  in  category  p,  one  may  calculate  a  coupling  matrix 
for  this  category.  The  coupling  matrix  is  related  to  the  directional  deriva¬ 
tive  of  the  probability  of  escape  P,(p,* d.  It. .), where  l. .  is  the  path 

•  Cijpij  1 J 

length  from  cell  i  to  cell  j  along  a  chord  inclined  to  the  radius  vector 
througja  cell  i  (Apruzese  ).  P^  is  a  nonlinear  function  of  the  path- 
integrated  optical  depth,  depending  upon  the  radiation  category's  damping 
mechanism  as  well  as  the  local  attenuation  parameter  values.  The  line 
radiation  categories  use  P^  functions  based  on  a  Doppler  line  profile,  while 
the  continuum  radiation  categories  use  a  P£  which  is  exponential  with  optical 
depth.  The  represent  geometric  projection  factors  needed  in  the  (dis¬ 
crete)  path  integral  estimate  and  can  be  calculated  for  all  radiation  cate¬ 
gories  once  the  mesh  *r^  is  given.  The  coupling  matrix  can  therefore  be 
defined  as 

+ 

X  ^  m 

C1jp  "  lVp’ f  *««(*)))  -  PE(P./  Jdxt(<i(x)))  , 

x(  X, 

where  x  is  a  path  distance  and  Ax^  *  x^  -Xj  is  the  cell  thickness,  along 
the  inclined  chord.  In  terms  of  the  coupling  matrix,  the  energy  lost  or 
gained  by  a  cell  is  the  inner  product 

Qphoto  j  *  ^  ^cepj  *  ?  Cijp  c  epi ^  ; 

while  the  energy  lost  from  the  plasma,  the  observable  output  from  any  cell 
in  any  radiation  category,  is  the  residual  sum 

Qout,jp  *  ce  pj  *  (1  '  ^  C2jp}  * 
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These  operations  are  done  with  the  aid  of  three  subroutines, 

SPFACX,  SPHERX,  and  PINT  due  to  R.  W.  Clark^  .  The  first  computes  the 

projection  factors  from  the  given  mesh  *r.  for  any  particular  angle  of 

•  J  J 

inclination  of  the  chord  through  each  cell.  The  second  and  third  act  in 

concert  to  produce  the  coupling  matrix  from  the  data  *d.  and  £...  The 
•  c  jp  1 j 

subroutine  PINT  produces  required  path  Integrals  of  escape  probability  in 

either  Doppler  or  exponential  mode,  while  SPHERX  produces  the  coupling  matrix. 

The  subroutine  TDOT  contains  all  of  these  calculations  but  provides 

options  for  bypassing  those  not  needed  for  the  iterative  refinement  of  the 

magnetic  field.  These  mode  switches  are  contained  in  the  common  block 

TRANSPORT-CHOICE  (along  with  the  CONTROL  matrix  for  subroutine  BRAGINSKII) 

as  the  logical  variables:  PHOTONS,  JSET-ONLY ,  EOSMODE,  and  DBLE-DOT-MODE. 

Of  this  set,  DBLE-DOT-MODE  is  referenced  by  CREOS  as  well.  The  flow  of 

the  subroutine  TDOT  is  illustrated  in  Fig.  5,  and  is  essentially  the  same 

sequence  of  operations  discussed  above  if  no  bypasses  are  taken.  The 

reared  Input  is  t*9ej,  JT„.  J  r, ,  *n,  j  ,  Jr, ,  *rjt  JEj ,  *Bj}. 

The  output  is  divided  among  the  common  blocks  TDOTOUT,  containing: 

{tn.  . ,  *t  j,  *9  .,  and  *0,,  ,  .  >  and  GRADS-FIELDS ,  containing 

c  i,o  c  i,j  c  e,j  c  e,j  rad,jp 

temperature  gradients,  compression  rates,  thermoelectric  fields,  total  plasma 
electric  fields,  classical  or  limited  axial  current  densities,  and  (if 
DBLE-DOT-MODE  is  .TRUE.)  information  required  by  MESHTRESS  and  HERTZDOTS 
for  higher  order  time  derivatives. 


FIGURE  5  -  Flow  of  Subroutine  TCOT 


6.  The  Subcycle  Calculating  Thermal  and  Electrodiffusive  Evolution 
in  the  Model  Plasma 

A  subcycle  is  preferred  for  thermal  conduction  (and  other  dissi¬ 
pative  processes)  because  the  parabolic  p.d.e. describing  the  evolution  is 
often  treated  best  by  an  implicit  method,  while  the  momentum  advance  can 
often  be  adequately  treated  by  an  explicit  method.  Also,  from  a  physical 
viewpoint,  the  fastest  timescales  of  the  model  plasma  are  those  associated 
with  electron/ion  thermal  exchange,  (unmagnetized)  electron  thermal  conduction 
and  radiative  losses  (under  conditions  of  strong  compression).  The  solution 
of  the  thermal  evolution  relations  (I II. 4.  a,  b)  for  a  fixed  level  of 
accuracy  will  therefore  require  more  time  derivative  evaluations  than  an 
equivalent  solution  of  the  velocity  field  relation  (III. 3),  even  allowing 
for  some  iterative  refinement  of  the  mesh  motion.  Finally,  as  a  matter  of 
flexibility  In  algorithm  choice,  the  subcycle  allows  one  to  separate  the 
thermal  evolution  from  the  mesh  evolution  within  the  code  and,  therefore,  to 
select  possibly  distinct  time-advance  mechanisms  appropriate  to  each. 

The  thermal  subcycle  presently  implemented  derives  from  the  observa¬ 
tion  that  a  spacetime  .  p.d.e.  which  has  the  derivative  operations  represented 
by  discrete  differences  on  one  domain  of  dependence  appears  as  a  coupled 
set  of  o.d.e.  on  the  orthogonal  domain.  The  particular  case  here  is  that  of 
spatial  differential  operators,  corresponding  to  finite  difference  operators 
(derived  through  conventional  or  smooth  interpolants) ,  producing  a  set  of 
coupled  equations  on  the  time  domain.  In  all  such  cases  the  coupling  of  the 
time  derivatives  Is  expressed  as  a  Jacobian  matrix,  implicitly  dependent  on 

the  discrete  mesh  used  for  the  derivative  operators.  Here  the  time  deriva- 

•  •  • 

tlves  are  the  material  derivatives  (T  ,  8  ,  E  }  and  the  underlying  mesh  motion 

C  C  4 

(r,  f,  F,  r*}  is  completely  transparent  to  the  subcycling  algorithm.  The 


three  fundamental  codes  for  this  process  are  described  below. 


(i)  GEARBOX,  the  subcycle  integrator 

Based  on  a  multi -time-level  variable  order  scheme,  the  subroutine 
GEARBOX  performs  the  thermal  subcycle  on  a  time  domain  [Tlo,  Thi]  defined 
externally.  The  input  consists  of  *Y.  *  V.  .,  2.0  . ,  VT  .... 

l0e  data  for  a  total  of  2*J  coupled  equations  of  motion. 

The  time-dependent  mesh  (with  which  the  thermal  fluxes,  congressional  heating 
and  radiative  cooling  are  computed)  is  specified  implicitly  over  this  interval 
for  Tlo  <  t*  <  Thi.  The  forward  integration  of  (Y.}is  done  by  the  IMSL 

J 

subroutine  OGEAR,  using  the  material  derivative  subroutine  FLUIDOTS  and  the 
coupling  matrix  (Jacobian)  produced  by  JACOB.  In  the  environment  seen  by 
DGEAR.the  problem  is  specified  completely  by  the  input  { 1 Y . > ,  the  interior 

J 

material  derivatives  {*'?.},  the  coupling  matrix  *(3Y -/3Y . ) ,  and  the  inter- 

J  *  J 

mediate  time  points  t*.  Values  of  t*  are  chosen  by  the  integrator,  in 
accordance  with  the  rates  of  change  produced  by  FLUIDOTS,  and  these  values 
are  thus  the  sole  input  arguments  available  for  the  specification  of  the 
Lagrangian  mesh  *r.. 

J 

To  sumnarize  the  operation,  GEARBOX  begins  by  checking  for  an  upper 
limit  in  the  explicit  times tep  supplied  as  input  whenever  the  integrator  is 
being  initialized.  Next  the  explicit  time  limits  set  as  arguments  are 
stored  in  the  common  block  <TIMEBASE>  for  use  by  FLUIDOTS  and  JACOB  in 
implicitly  advancing  the  fluid  mesh.  If  desired,  a  short  report  on  the  model 
thermodynamic  variables  is  written.  The  actual  Integration  is  effected  by 
a  call  to  OGEAR,  and  upon  completion  various  error  flags  are  checked  and  the 
severity  of  any  integrator  errors  is  assessed.  If  desired  another  short 
report  on  the  evaluation  is  produced;  and,  if  a  terminal  integrator  error 
Is  found  a  completed  dump  of  all  pertinent  Information  Is  triggered.  The 
subroutine  returns  with  a  new  (provisional )  vector  {i+1Y.}  overwriting 
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the  original  Input.  The  advanced  temperatures  are  either  accepted  or 
iterated  upon  in  the  calling  routine,  depending  on  the  chosen  option  for  the 
fluid  mesh  evolution. 

The  integrator  DGEAR  offers  several  control  parameters  and  con¬ 
siderable  flexibility  in  the  integration  scheme.  The  most  important  options 
are  the  use  of  explicit  Jacobians  or  estimated  Jacobians,  and  the  control 
of  the  local  truncation  error  through  a  variable  TOL.  The  explicit  Jacobian 
is  apparently  the  most  favorable  option  in  this  hydrodynamic  application,  but 
if  the  application  of  more  detailed  physics  demands  a  very  expensive  evalua¬ 
tion  for  this  matrix,  an  estimate  based  only  on  finite  differences  is  adequate. 

The  variable  TOL  provides  a  natural  and  consistent  means  of  specifying  time 
integration  parameters  for  the  mesh  evolution  as  well,  cf.  Sections  (.III .H . i  &  ii) 
The  integration  scheme  can  be  selected  as  either  an  implicit  Adams  method  of 
up  to  twelfth  order  or  a  backward  differentiation  method  of  up  to  fifth  order 
(Gear's  stiff  method).  Both  methods  are  of  the  implicit  linear  multilevel 
type  and  require  the  solution  of  an  algebraic  system  at  each  interior  (sub¬ 
cycle)  timestep.  The  basis  for  an  optimal  choice  between  the  two  lies  in 
the  performance  of  several  benchmark  calculations  discussed  in  Section  III.  I 
below. 

(ii)  FLUIDOTS,  the  material  derivatives 

The  physical  content  of  the  problem  is  defined  entirely  in  the 
subroutine  FLUIDOTS,  which  incorporates  an  explicit  time  advance  of  the 
data  (*Zfe,  *3TZkf*rj}  in  order  to  define  the  mesh  and  electromagnetic  fields 
required  to  compute  In  summary,  FLUIDOTS  receives  as  input  (*Yj)  and 

t*,  which  is  any  Interior  point  selected  by  DGEAR  on  the  <TIMEBASE>  for  the 
calculation  of  a  derivative.  First  FLUIDOTS  maps  *Y.  to  separate  variables 

w 

*T"lj’  Next  the  mesh  is  obtained  by  a  forward  Taylor  series  using 


whatever  Information  resides  In  the  common  block  MESHADVDATA: 


*r.  =  V.  +  (t*  -  Tlo)1^  +  -5-  (r*  -  Tlo)2  V.  +  Ut*  -  Tlo)3  %r. 
J  J  j  <■  Jo  J 

*r.  -  +  (t*  -  Tlo)1^  +  £{t*  -  Tlo)2  Vj  . 


This  mesh  Is  checked  for  inversion  {*r^.  >  *pj*i}  and  if  such  inversion 
occurs  a  terminal  error  condition  is  generated  and  a  dump  of  all  pertinent 


information  is  obtained.  The  mesh  is  also  checked  for  closure  of  the  annulus 


(*r,  £0)  and,  if  closure  is  sensed  during  a  subcycle,  the  motion  of  the 
innermost  cell  boundary  is  corrected,  *r^  -*»  0  &  *fj  *  0.  Once  these  opera¬ 
tions  are  complete  an  array  of  cell  centers  is  generated  by  equal  areal 
partition  and  a  central  velocity  *r.  is  assigned  by  area  weighting.  A  cell 

C  j 

ion  density  is  also  calculated  from  the  (conserved)  particles  resident  in 
the  cell. 


The  switches  EOSMODE  and  JSET-ONLY  are  then  set  for  a  complete 


TOOT  calculation.  If  electromagnetic  effects  are  turned  off  this  single 

TDOT  calculation  proceeds.  If,  as  is  usually  the  case,  the  E,  and  B.  are 

needed  then  JSET-ONLY  is  set  .TRUE,  and  ZPUSH  advances  the  Hertz  potential 

from  its  last  data  base  *  *3  Z.)  to  its  values  at  t*.  Given  the  results 

of  ZPUSH,  the  subroutines  EFIELD  and  BFIELD  are  called  to  establish  Ez 

completely  and  partially,  using  the  current  densities  of  the  previous 

x*  point.  Using  the  provisional  *Bflj,  a  call  to  TDOT  with  EOSMODE  and 

c  wj 

JSET-ONLY  set  .TRUE,  recalculates  the  central  current  densities  at  fixed 


Ez  and  estimated  Z.  A  second  call  to  BFIELD  uses  the  new  current  densities 
and  thus  implicitly  refines  1(B).  After  one  refinement  pass  (TDOT; 

BFIELD),  EOSMODE  is  set  .FALSE,  and  the  initial  calculation  of  j , 
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and  *Tgj  remains  in  place.  The  gradient  calculation  remains  undisturbed 
as  well  so  that  the  Ez  thermoelectric  field  component  will  change  only 
through  the  refinement  of  its  BQ-dependent  coefficient.  The  refinement 
passes  continue  until  {£8Qj}  convergence  is  achieved  to  a  specified  tolerance 
or  until  a  maximum  number  of  passes  has  been  made.  Then  JSET-ONLY  is  set 
.FALSE,  and  a  final  call  to  TOOT  calculates  the  full  complement  of  material 
derivatives  and  transport  coefficients  using  the  self-consistent  magnetic 
field  values.  In  this  final  pass  the  gradients  are  recomputed  because  the 
thermal  conduction  has  changed  with  the  magnetic  field  refinement.  If 

magnetic  field  convergence  is  not  ach levee  to  the  desired  accuracy,  e.g. 

-4 

o(10  ),  the  vector  of  fractional  changes  at  the  last  iteration  is  output  as 
a  diagnostic,  and  if  desired  a  full  report  of  the  calculation  is  output 
as  well. 

The  last  operations  performed  by  FLUIDOTS  are  to  map  *9  . ,  *tT  . 
to  the  single  variable  and  to  update  the  acceleration  and  jerk  arrays 

,  Vj  If  required.  An  update  of  these  variables,  contained  in <MESHADVDATA>, 
Is  triggered  within  the  subcycle  whenever  t*  is  larger  than  any  previous 
time  argument  requested  T*>rlast,  or  If  t*  ■  Thl.  As  discussed  in  the  following 
Section  III.H,  this  enables  one  either  to  correct  the  mesh  evolution 
smoothly  as  temperatures  evolve  or  to  set  up  for  a  complete  iterative 
refinement  of  the  mesh  evolution.  The  mesh  update  occurs  with  a  call  to 
MESHTRESS  and.  If  OBLE-OOT-MOOE  Is  .FALSE,  as  is  the  case  during  a  sub¬ 
cycle,  It  Involves  only  the  explicit  recalculation  of  the  acceleration 
array  *Pj.  The  jerk  array  is  continually  updated  as 


for  t*  >  T-|ast.  When  FLUIDOTS  is  called  outside  a  subcycle,  DBIE-OOT-MOOE 


is  set  .TRUE,  and  an  explicit  evaluation  of  the  jerk  array  Vj  is  obtained 
as  well. 


(iii)  JACOB,  the  link  to  the  original  p.d.e. 

The  coupling  matrix  calculated  by  JACOB  is  a  simple  if  lengthy 
exercise  in  the  use  of  partial  derivative  chain  rules.  In  the  case  of  a 
simple  three  point  difference  formulae  this  matrix  is  slightly  sparse  and 
pentadiagonal .  The  compression  and  exchange  heating  create  a  block-tri diagonal 
band  dependent  on  (V  • V )  and  the  local  Tg,  while  the  thermal  conduction 
operators  add  two  more  bands  and  some  diagonal  contributions.  In  the  case 
of  a  smooth  interpolant  (G  or  E-spline)  the  matrix  contains  more  nontrivial 
bands  of  monotonically  decreasing  importance  away  from  the  diagonal.  The 
rate  of  decay  is  determined  by  the  coupling  parameter  choice  in  the  under¬ 
lying  interpolant. 

However  one  defines  the  differentiation  process,  the  general  form 
of  the  Jacobian  is  easily  derived  from  the  structure  of  the  subcycled  array 
{Yj}.  For  all  entries  involving  the  material  derivative  of  9fi  the  local 
value  of  the  branching  ratio  is  the  element  of  central  interest,  i.e. 


3  In  b\ 


9 


while  only  on  the  first  subdiagonal  does  0  couple  to  Tj,  since 


apart  from  some  radiative  contributions  due  to  the  cell's  self-opacity. 
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In  a  similar  manner  the  entries  involving  the  material  derivative  of  Tj 
couple  to  8  only  on  the  first  super-diagonal,  i.e. 


III.  a  (]  .  Ua  »)  .  Ts- 

seJ  bj  A  Ve.j 


The  coupling  on  the  main  diagonal  consists  of  several  contributions.  The 
primary  one  is  congressional  heating  (-  |  V‘V),  with  exchange  losses  and  the 
term  from  the  central  coefficient  of  whatever  gradient  formula  is  employed 
in  computing  the  heat  flux  also  important.  In  the  case  of  9g  a  radiative 
diagonal  contribution  also  arises. 

The  general  Jacobian  structure  obtained  for  a  three-point 
difference  scheme  is  illustrated  in  Fig.  6.  There  compressional  heating  is 
denoted  by  9  on  the  diagonal;  the  ion-electron  exchange  terms  are  shown  by 
x;  and  the  contributions  from  the  spatial  differencing  in  thermal  con¬ 
duction  are  denoted  by  A. 


Figure  6  -  A  simple  Jacobian  matrix 
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The  inclusion  of  radiation  terms  is  dependent  on  the  coupling 
matrix  simulations  discussed  above  in  computing  Qphoto-  It  is  a  straight¬ 
forward  but  more  complex  calculation  along  these  same  lines. 


Civ)  Modifications  for  Ez  diffusion 

If  the  option  of  diffusive  evolution  for  Ez  is  desired,  an 
alternate  trio  of  subroutines  is  employed.  The  vector  Y.  now  becomes 

J 

{.  *Tt4,  *E  .,  in  analogy  to  the  structure  used  in  GEARBOX. 

This  implies  a  Jacobian  which  is  block  pentadiagonal  with  further  bands 
added  depending  on  the  spatial  differencing.  The  routine  GEARBOX  is 
replaced  by  TETGEAR;  FlUIDOTS,  by  TETDOTS;  and  JACOB,  by  TETJAC.  The 
routine  TETDOTS  contains  most  of  the  principal  differences;  it  must  call 
an  additional  subroutine  EDOT  which  evaluates  equation  (III. 2).  The 
Iterative  refinement  of  B0  appearing  in  FLUIDOTS  persists,  but  the 
external  circuit  equation  requires  an  iterative  refinement  of  and 
c^zj  together  in  order  to  achieve  a  self-consistent  prescription  of  the 
boundary  conditions  for  equation  (III. 2).  Once  the  values  of  *f  , 

are  established  by  TDOT  and  MESHTRESS,  however,  all  iterative  refine¬ 
ment  can  be  made  internal  to  EDOT.  The  flow  of  TETDOTS  is  thus  quite  similar 
to  that  of  FLUIDOTS  with  an  additional  call  to  EDOT  prior  to  the  final 
return. 


H.  The  Complete  Fluid  Evolution  Package  -  HYDROPUSH 

All  of  the  foregoing  algorithms  are  combined  in  the  subroutine 
HYDROPUSH  .which  forms  the  nucleus  of  Program  ZDIPR-apart  from  startup, 
diagnostic  and  graphics  modules.  HYDROPUSH  is  concerned  exclusively  with  the 


advance  of  the  <FUJID-$TATE>  over  the  variable  major  time  step  which  it 
selects  as  appropriate.  The  only  required  inputs  are  the  common  blocks 
<FLUID-STATE>  and  <FIELDADVDATA>  (containing  the  previous  time  slice  of 
fluid  variables  and  field  variables)  and  the  appropriate  mesh  dimension, 
scale  factors  and  physical  constants, obtained  from  a  variety  of  sources. 

The  output  is  an  update  of  <FLUIDSTATE> ,  and  an  increment  of  the  variable 
MSI  (Main  Step  Index)  by  1.  If  a  restart  file  is  desired,  it  is  created 
for  later  use  and  assigned  a  record  number  equal  to  the  MSI.  If  any  severe 
errors  occur  in  the  advance,  dump  files  are  created  to  allow  examination 
of  intermediate  results,  and  a  variety  of  reports  at  intermediate  phases 
of  the  calculation  are  available  as  general  diagnostics.  The  architecture 
of  HYDROPUSH  assumes  that,  after  the  fluid  advance  is  successful,  the 
main  program  will  call  ZPUSH  separately  to  advance  the  fields  1 Z k  and 
^3TZk  -  unless  the  electrodiffusive  approximation  is  being  used.  In  that 
case  HYDROPUSH  will  also  advance  the  E.  field  by  means  of  TETGEAR. 

The  sequence  of  processing  begins  by  reading  elements  of  the 
<FLUID-STATE>  into  <MESHAOVDATA>  and  the  subcycle  vector  { 1' Y . } .  Once 
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the  appropriate  data  base  is  inferred  a  full  complement  of  material  deriva¬ 
tives  Is  computed  by  FLUIDOTS/HERTZDOTS  or  TETDOTS.  These  material  deriva¬ 
tives,  including  the  field  ’(all  J  ).  or  ].£  .  ,  are  then  used  by  STEPPER  to 
select  a  major  time  step.  The  environment  of  either  GEARBOX  or  TETGEAR 
is  held  if  one  wishes  to  iterate  the  mesh  evolution  (GRIDITER=.TRUE. ) . 
Otherwise  a  single  subcycle  is  done  (advancing  1+1Y.,  i+1r.  and  1+1r.) 

J  J  J 

and  the  new  FLUID-STATE  is  output. 

The  iterative  refinement  of  the  mesh  evolution  is  accomplished  by 
first  making  a  predictor  step  with  the  jerk  continuously  updated  over  the 


subcycle.  At  the  end  of  the  first  pass  the  GEARBOX  or  TETGEAR  environ¬ 
ment  is  reset  to  its  content  at  the  beginning  of  the  subcycle.  Then  the 
full  compliment  of  material  derivatives  is  explicitly  extracted  from  the 
(provisional)  advanced  fluid  data,  and  the  jerk  coefficient  is  overwritten 
by  a  linear  combination  of  its  values  at  the  beginning  and  end  of  the  step. 
Further  passes  over  the  subcycle  use  progressively  better  estimates  of  the 
jerk  until  convergence  in  the  linear  combination  coefficients  is  achieved 
or  the  maximum  allowed  number  of  corrective  passes  is  reached.  A  final 
pass  with  the  fully  self  consistent  MESHADVDATA  is  then  executed  to  obtain 
the  most  accurate  thermal  evolution.  The  resulting,  time-advanced  GEARBOX/ 
TETGEAR  environment  is  left  in  place  and  the  new  <FLUID-STATE>  is  output. 

The  general  structure  of  HYOROPUSH  is  illustrated  in  Fig.  7;  the 
details  of  the  time  step  algorithms  and  the  iterative  refinement  scheme 
are  discussed  In  the  subsections  below. 

(i)  Time  step  considerations 

Over  a  subcycle  the  steps  are  controlled  by  DGEAR,  but  the  assign¬ 
ment  of  the  proper  <TIMEBASE>  must  be  based  on  those  rates  not  accessible 
to  the  subcycle  Integrator,  viz.  {^r.,,  V.*,  V.,  *(3"j  ).  or  *6  >.  The 
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comparison  among  these  rates  ,  and  those  time  steps  derivable  from  them,  is 
accomplished  by  the  subroutine  STEPPER. 

The  velocities,  accelerations,  and  jerks  of  the  mesh  can  be  used 
to  formulate  several  distinct  time  step  estimates.  First  is  the  magneto- 
sonic  Courant-Fredrich-lewy  (CFL)  condition,  which  seeks  to  insure  that  the 
spatial  domain  of  influence  on  future  values  is  contained  in  the  domain 
of  dependence  established  by  the  finite  differencing  algorithms  which  pro¬ 
vide  the  acceleration.  The  domain  of  influence  is  specified  by  the  most 
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U|Ml4t«  <HUin  STAH> 


rapid  signal  velocity  physically  significant  in  the  fluid,  which  for  the 
present  application  is  the  magnetosound  speed 


The  domain  of  dependence  for  the  pressure  gradient  stresses  is  at  least  two 
complete  cells,  so  the  distance  of  propagation  is  one  cell  width.  Since  the 
cell  is  generally  either  compressing  or  expanding,  the  disturbance  propaga¬ 
tion  is  superposed  on  the  relative  boundary  motion,  and 


!™  *5i(Vjn  -  S'/'c’Vj 


with  6  i  in  the  range  j -►  j  for  explicit  mesh  integration  and  perhaps  as 
2  7 

large  as  j  -*•  ^  for  implicit  evolution. 

A  second  criterion  is  the  Local  Truncation  Condition  (LTC)  based 
upon  the  highest  time  derivative  of  the  forward  Taylor  series  used  in 
advancing  the  mesh  f*j.  Since  area  weighed  differencing  is  the  basis  of  the 
gradient  calculations  the  highest  time  derivative’s  contribution  to  the 
relative  cell  area  change  should  be  (approximately)  the  truncation  error 
^GEAR  the  inte9rator  (the  GEARBOX  variable  TOL)  multiplied  by  the 
typical  number  of  subcycles  expected,  i.e. 

4tj  ' LTC  '  jv(  “j+rV'  j  ^ 

with«2A.  (5-io)  •  aGEAR. 

A  third  related  criterion  Is  the  relative  cell  area  change  itself, 
calculated  through  second  order  in  the  Taylor  series  -  the  Quadratic  Cell 
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Area  change  (QCA).  With  AOj  =  and  similar  notation  for  the 

material  derivatives,  an  expanding  cell  will  grow  by  a  fraction  £  in 


i  -A 

QCA  ^ 


/ 2  sW. 

n  +  -3 - i — 1  -  i 


where  6  may  be  set  in  the  range  0.08  -*■  0.20.  A  similar  expression  exists 

for  the  fractional  shrinkage  of  a  contracting  cell,  and  the^At  J  }  is 

J  QCA 

just  the  collection  of  such  estimates.  A  cell  destined  to  change  from 
contraction  to  expansion  (or  vice  versa)  over  the  time  step  is  excluded 
from  the  selection  algorithm. 

A  fourth  mesh-based  criterion  is  the  time  for  collapse  on  axis. 
If  the  innermost  zones  (r^  and  r^)  are  projected  ahead  using  the  Taylor 
series  coefficients,  then  one  seeks  a  time  which  allows  the  first  to 
intercept  the  axis,  but  not  the  second.  After  the  first  has  collapsed, 
only  the  time  for  axial  interception  by  th*  second  is  relevant.  The 
collapse  time  ^At  Is  then  appropriately  set  in  between  these  two. 

v* 

Finally  a  step  £t|  can  be  set  by  insisting  that  the  relative 

£JH 

changes  in  wave  sources  or  diffusion  sources,  i.e.  3tJ,  E  be  bounded  by  a 
fraction  similar  to  that  limiting  cell  area  changes. 

The  subroutine  STEPPER  calculates  these  various  time  steps 
{At|CpL,At|ltc,At|qCA,At|c  ,  At  >  and  then  assigns  a  final  At  as 


ilinf  ^At.. 
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,  inf  1 At^ |  ,  inf  nAT 

CFL  j  J  LTC  j 


,  |  ,  At  , 
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This  provisional  At  is  then  checked  for  any  danger  of  mesh  inversion  by 
advancing  the  mesh  to  the  end  of  the  interval.  If  a  mesh  inversion  is 
implied  the  step  size  is  cut  by  a  factor  of  0.9  and  the  test  is  made 
again.  This  procedure  continues  until  either  the  step  cuts  prevent  the 
Inversion  or  the  original  step  has  been  cut  to  a  tenth  of  its  original 
size.  In  the  later  event  the  calculation  proceeds  with  a  warning  of  possible 
inversion.  STEPPER  then  outputs  the  new  end  point  (Thi)  for  the  subcycle 
integrator  and  an  estimate  (delTmin)  to  be  used  as  an  explicit  step  if  the 
time  stack  internal  to  the  DGEAR  integrator  is  to  be  reinitialized. 


(ii)  Explicit  and  implicit  mesh  advance  with  HYDROPUSH 

The  simplest  mesh  evolution  available  is  an  explicit  advance  using 
the  Taylor  series  coefficients  at  a  single  time  slice, 


1 


„  +  *A  ii  +  *At2  i..  .  *At; 

r j  +  At  j  ~T~  rj  +  — 


i. 


Fj 


ir .  +  *AT1f.  +  V. 


The  local  truncation  condition  in  STEPPER  insures  that  this  Taylor  series 
receives  a  small  contribution  from  the  highest  order  (and  potentially 
most  noisy)  term  V^.  In  order  to  minimize  the  global  accretion  of  what¬ 
ever  error  might  be  added  in  any  one  such  time  step,  however,  it  is 
necessary  to  insist  on  an  evolutionary  self-consistency  in  the  advance  of 
the  acceleration  fields  and  the  jerk  fields.  4  consistency  determined  at  the 
future  <FIUID-STATE>  by  the  position,  velocity  and  temperature  fields  and 
external  stresses. 
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One  way  to  help  provide  self-consistency  is  to  alter  r*.  during 

J 

the  subcycle  by  examining  the  acceleration  values  obtained  (*r.)  and 

<1 

revising  by  the  rule, 

<-  *r,  -  V , 

★y*  3  J  -1 

j  *At 


continuously  as  the  subcycle  progresses,  the  superscript  t-  denotes  a 
forward  time  ordering.  This  is  tantamount  to  choosing  a  (variable  interval) 
time  averaged  for  the  mesh  advance,  and  clearly  has  a  cumulative  but 
distinct  impact  over  the  subcycle  on  the  position  and  velocity  fields. 

For  the  grid  *r.  this  version  Implies  a  local  effective  acceleration 

J 


^  j*  +  j.  (*i* 

ri  3  1  rj 


9 


while  for  the  velocity,  the  local  effective  acceleration  is 


These  effective  accelerations  appear  to  be  available  only  for  future 
subcycle  times  -  *t*j  is  the  acceleration  derived  from  the  fluid  upon 
arri val  at  t*  by  means  of  the  Taylor  series  coefficients  previously 
computed.  The  GEARBOX  integrator  is  not  always  committed  to  a  monotone 
increasing  march  in  t*  however,  so  that  any  mesh  required  when  it  "backs 
up"  should  reflect  the  forwardmost  Taylor  series  information  available. 
For  this  reason  the  material  derivative  codes  FLUIDOTS  and  TETDOTS  keep 
track  of  the  times  t*  input  to  them  and  perform  the  revision  of  T.  (and 

J 

hence  the  calculation  of  the  effective  accelerations)  if  and  only  if 
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the  given  t*  is  the  largest  time  argument  seen  thus  far  in  the  subcycle. 

This  approach  therefore  continuously  updates  the  Taylor  series  on 
an  ever  broader  temporal  data  base  in  order  to  follow  the  general  trend  of 
Tj  instead  of  just  a  single  value.  The  explicit  (or  predictor)  phase  of 
HYDROPUSH  is  just  the  implementation  of  this  algorithm. 

On  the  other  hand,  nothing  guarantees  that  the  smooth  changing  of 

Tj  is  sufficient  to  provide  true  self-consistency  in  the  motion.  At  the  end 

of  a  subcycle  there  exist  two  versions  of  The  first  is  available 

from  a  time  slice  of  advanced  fluid  variables  ,  i+1rj,  ^+cTIj’ 

,  ^+*EZ),  the  second  is  that  derived  over  the  time  step  using  the  jerk. 
-i i+lu  _  i  •  •  .  f*_  ^ 


coefficient 


rj- 


The  mean  value  theorem  guarantees  however  that  at  least  one  value 
of  Tj  (on  Tlo  <r*<_  Thi)  exists  which  makes  the  second  version  of  1+1fj 
exact  to  all  orders.  The  essential  thrust  of  the  implicit  method  described 
here  Is  to  find  that  best  value  for  each  cell  which  provides  a  solution 
to  the  nonlinear  difference  equation 


1+ltt  /1+1  1+1*  i+lT  i  +  lc  i+la  i+lp 

*jF(  ry  rj’  c'lj’  c9ej  ’  c  Ezj’  c  Bej  ’  cErj' 


'*SF  *  \ 


to  some  specified  precision,  e.g.  the  TOL  in  GEARBOX.  Here  the  accelera¬ 
tion  at  the  new  time  level  depends  on  the  entire  time  history  of  the  fluid 
variables  computed  over  the  subcycle.  The  most  general  method  of  solution 
is  to  iterate  across  the  subcycles.  This  is  accomplished  by  representing 
the  unknown  best  value as  a  weighted  sum  of  the  (explicitly  calculable) 
end  point  values,  viz. 


and  iteratively  computing  by  several  passes  across  the  subcycle  until 

convergence  (in  <5.)  is  achieved.  At  any  iteration  (p)  the  new  5.  is 
J  J 

computed  from  the  implied  new  fluid  state  (using  the  last  <5^  vector)  and  the 
original  fluid  state  by  the  required  solution  of  the  difference  equation,!* .e. 


psj 


(p-i  -  ^if)  . 


V 


At 


jF 


WL, 

p-l?jF 


-  V 


jF 


The  jerk  coefficient  (embedded  in  MESHADVDATA)  is  then  updated  to  be  the 
current  estimate  of  A‘fjp  given  above  and  the  entire  subcycle  is  recalculated, 
completing  the  loop.  Convergence  to  a  particular  vector  <5.  is  thus  equiva- 

J 

lent  to  a  solution  of  the  difference  equation  and,  in  the  limit  of  small  time 
steps,  any  solution  of  the  difference  equation  closely  approximates  a 
solution  of  the  original  partial  differential  equations.  In  the  implicit 
scheme  of  course  the  revision  of  T.  during  the  subcycle  is  suppressed. 

J 

The  present  implementation  of  this  scheme  in  HYDROPUSH  relies 
on  the  simple  replacement  algorithm  discussed  above  in  order  to  find  6.. 

J 

It  Is  probable  that  convergence  can  be  improved  by  using  more  sophisti¬ 
cated  estimates  based  on  the  iterative  sequence  {  5^}  p*0, 1,2,3..,  but  the 
best  such  choice  will  come  from  examining  the  algorithm's  performance 
and  thus  cannot  be  selected  a  priori  with  any  certainty.  It  is  also 
probable  that  the  Inclusion  of  previous  time  levels  (1-1,  i-2...)  in  an 
estimate  of  would  be  a  means  of  filtering  the  noise  that  might  develop 


In  this  variable,  and  general  methods  of  effecting  this  extension  are 
under  study. 


I.  The  Planned  Sequential  Benchmark  of  This  Method 

a)  Checks  on  internal  energy  conservation 

On  a  stationary  mesh  (no  congressional  heating  or  expansion 
cooling)  with  radiative  losses  disabled,  an  initial  separation  of  Tj  and 
Tg  (©e)  will  relax  inhomogeneous ly  when  the  (.time  invariant)  ion  density 
varies  in  space.  As  a  test  of  the  subcycling  algorithm.the  time  asymptotic 
final  state  (Tj  *  Tg  (0^)  »  for  all  values  of  r)  must  contain  the 
same  internal  energy  as  the  original.  Preliminary  tests  with  uniform  ion 
density  (and  uniform  relaxation)  have  demonstrated  this  energy  conserva¬ 
tion  to  0(10-3). 

b)  Checks  on  the  conservation  of  kinetic  and  internal  energy 

There  exists  a  class  of  self-similar  Gaussian  implosions  and 

explosions  for  which  the  time-dependent  density  scale  length  (or  any  other 
Lagrangian  marker  radius)  varies  as 

with  Yg  a  initial  thermal  energy/initial  flow  kinetic  energy.  These  flows 
remain  Isothermal  in  space  at  all  times  (removing  the  effects  of  thermal 
conduction)  and  heat/cool  8g  and  Tj  equally  (removing  the  ion/electron 
thermal  exchange.  It  Is  thus  a  necessary  check  to  insist  that  the  algorithm 
preserve  this  flow  locally  on  the  mesh  and  also  conserve  the  sum  of  flow 
and  internal  energy  globally  in  time. 

Once  these  benchmarks  are  obtained  the  numerical  performance  of  the 
code  can  be  properly  optimized  with  respect  to  mesh  density,  time  step,  and 
convergence  criteria. 
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c)  Checks  on  the  conservation  of  kinetic,  internal  and  magnetic 

field  energy 

For  a  strict  MHO  model  plasma  with  a  fixed  (constant)  current,  the 
self  similar  solutions  due  to  Felber16  provide  oscillatory  trajectory 
families  for  Gaussian  implosions.  One  makes  contact  with  these  in  the 
present  model  by  legislating  a  particular  current  profile,  and  by  neglecting 
ohmic  heating,  ambipolar  £r  stresses,  Ionization  dynamics  and  radiation. 

The  resulting  motion  can  then  be  compared  to  the  analytical  results  over 
several  oscillations,  and  the  energy  interchange  among  these  three  reservoirs 
monitored  with  regard  to  its  global  conservation  in  time. 

d)  Examining  the  process  of  collection 

A  final  basic  study  involves  the  filling  of  the  central  cavity  by 

an  initially  annular  load.  In  cylindrical  geometry  this  process  is  quite 

similar  to  the  diaphragm  problem  studied  by  A.  Lapidus17  and  a  detailed 

comparison  of  numerical  results  will  be  useful  in  assessing  the  shock  and 

rarefaction  resolution  capabilities  of  HYDROPUSH.  In  the  limit  of  large 

radii  the  annular  diaphragm  evolution  is  very  close  to  the  corresponding 

% 

planar  problem  for  which  analytic  flow  solutions  are  known. 


66 


References 


1.  R.  E.  Terry  and  J.  Guillory,  Development  and  Exploration  of  the  Core/ 
Corona  Model  of  Imploding  Plasma  Loads,  JAYCOR  report  number  TPD-200- 
80-001- DFR,  1980. 

2.  D.  G.  Colombant,  M.  Lampe  and  H.  W.  Bloomberg  WHYRAC,  A  New  Modular 
One-Dimensional  Exploding  Wire  Code,  NRL  Memo  Report  3726,  1978. 

3.  R.  E.  Terry  and  J.  Guillory,  Development  and  Exploration  of  the  Core/ 
Corona  Model  of  Imploding  Plasma  Loads,  JAYCOR  report  number  TPD-200- 
80-001- FR,  p.92. 

4.  ibid  p.  44  and  45,  pp  95-98. 

5.  R.  E.  Terry  and  J.  Guillory,  Annual  Report  on  Modeling  of  Imploding 
Annular  Loads,  JAYCOR  report  number  J207-81-004,  1981,  p.  20. 

6.  ibid,  pp  26-58. 

7.  Ibid,  Appendix  0. 

8.  R.  E.  Terry,  J.  Guillory  and  D.  Duston,  Generalized  Hertz  Vector 
Potentials  and  Their  Application  to  Diode  Imploded  Plasmas,  Bull.  APS, 
26,  No.  7,  Sept.  1981. 

9.  A  Talml  and  G.  Gllat,  Journal  of  Computational  Physics  23,  p.  93,  1977. 

10.  D.  Duston  and  J.  Davis,  Phys.  Rev.  A  21.,  May  1980. 

11.  Y.  T.  Lee  (private  communication) . 

12.  R.  E.  Terry  and  J.  Guillory,  Development  and  Exploration  of  the  Core/ 
Corona  Model  of  Imploding  Plasma  Loads,  JAYCOR  report  number  TPD-200-80 
001-FR,  1980,  p.  70. 

13.  S.  I.  Braglnskli  in  Reviews  of  Plasma  Physics,  ed.  M.  A.  Leontovich 
(Consultants  Bureau,  NY  1965). 


67 


J.  A.  Pruzese,  J.  Quant.  Spectrosc.  Radi  at.  Transfer  £5,  p.  419,  1981. 

R.  Clark  (private  communication). 

F.  S.  Felber,  Self-similar  Oscillations  of  a  z-pinch  (to  appear  in  Physics 
of  Fluids). 

A.  Lapidus,  J.  Computational  Phys.  8,  p.  106,  1971. 


APPENDIX  I!  VAX  INTERACTIVE  INTERPOLATION  PACKAGE 

A«  an  illustration  of  the  coeoands  REPRESENT  and  EVALUATE ,  the 
fcllouina  calculation  produced  tho  G-srline  coofficionts  of  tho  branch- 
infl  ratio  in  Aluoinum  cf .  Charter  II. 

Tho  fllo  of  input  data  reside*  on  <CHEMNODE>*  I.E. 


PILES  CHEHNODE 

•1  rectory  -0RA1 1CIPRI. HYDRO . CHEM3 


CRS1S.DAT 11 

7/11 

12-AUS-1981 

18142 

(RE* RUE, * 

CRE17.DAT 11 

7/11 

7-AUG-1981 

09130 

(RE, RUE,* 

CRE19 , DAT 1 1 

7/11 

19-AUG-1981 

13112 

(REtRUE,* 

CRE19DR.DAT 11 

S/ll 

20-AUG-1981 

14141 

(RE* RUE*  * 

CRE19LTED.DAT1 1 

7/11 

24-MAR- 19S2 

14148 

(RE, RUE, * 

CRE19LZTS.DAT11 

7/11 

20-AUG-1981 

14137 

(RE* RUE, * 

CRE19ZE.DAT! 1 

7/11 

20-AUG-1981 

14117 

(RE* RUE*  * 

CRE21.DAT 11 

7/11 

7-AU0-1981 

09132 

(RE* RUE, * 

CRECLAM.DAT 11 

4/11 

23-MAR- 1982 

20148 

(RE* RUE, * 

CREKLAM.DAT 11 

3/11 

19-MAR- 1982 

10132 

(RE* RUE, * 

CRELLAM.DAT 1 1 

4/11 

23-MAR- 1982 

13139 

( RE » RUE , , 

CRERLAM.DAT 1 1 

4/11 

24-MAR- 1982 

13124 

(RE  *  RUE , * 

EHITALC.DAT! 1 

3/11 

23-MAR- 1982 

20144 

(RE*RUE*  * 

EH I TALK. OAT I 1 

3/11 

19-MAR-1982 

10121 

(REtRUE** 

ENITALL.DAT 11 

3/11 

23-HAR-1982 

13129 

(REtRUE, , 

EHITFILE.EXE! 1 

19/22 

17-MAR-1982 

14119 

(RE  *  RUE , * 

E0SF1LE.EXE! 1 

8/11 

19-AU8- 1981 

11119 

(RE* RUE, * 

EOSFILCXP.EXE 11 

10/11 

19-AU0-1981 

11120 

(REtRUE,* 

RCR8ECT.DAT 11 

3/11 

17-MAR- 1982 

14120 

(REtRUE** 

Total  of  19  files*  121/220  blocks. 

and  is  naaod  CRE19LZTS.DAT 1 1 .  Tho  coaaand  REPRESENT  IS 

SHOW  BY  MOL  REPRESENT 
REPRESENT  ■  (COMMANDS ! REP IT • CON 

TYPE  COMMAND*: REP IT. COM 


««  REPRESENT  P1LENAME  »» 

PI 

ASSIGN  'PI '.DAT  INFILE 

INQUIRE  P2  *UHAT  NAME  DO  YOU  U1SH  FOR  THE  OUTPUT  FILE?  ENTER  IT, PLEASE.’ 
ASSIGN  'P2' . DAT  OUTFILE 
ASSIGN  ' P2 ' H «  DAT  HEADER 
ASSION/USER.HODE  SYSt COMMAND!  SYSt INPUT’. 

RUN  UTlLITYlREPtT.EXE 

!««  PROMPTS  UILL  BE  GENERATED  FOR  FURTHER  INFORMATION.  »» 

DEASSION  INFILE 
DEASSION  OUTFILE 
DEASSION  HEADER 
FILES  t.DAT 


and  its  use  is  indicated  by  tho  above  coooand  file. 
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REPRESENT  CHEMNODE ! CRE19LZTB 

UNAT  NAME  00  YOU  WISH  FOR  THE  OUTPUT  FILET  ENTER  IT, PLEASE. 5  LAHZTB1? 

UMAT  IS  YOUR  SPLINE  CHOICE? 

ENTER  'QSPLINES. ' , ‘ESPLINES. * ,0R  'EI2SPLINES. * 

'OSPLINESi  ' 

00  YOU  WISH  INSL  ACCURACY  CHECKS? 

'YES' 

ENTER  THE  DI6XTS  OF  PRECISION  FOR  YOUR  INPUT  FILE. 

S 

00  YOU  WISH  INSL  ITERATIVE  REFINEMENT? 

'NO  ' 

ARE  THERE  BOUNDARY  DERIVATIVES  IN  THE  INPUT  FILE? 

'NO' 

ENTER  THE  NUMBER  OF  X-VALUES  IN  YOUR  INPUTFILE. 

31 

ENTER  THE  NUMBER  OF  FUNCTIONS  (  F<X>>  TO  BE  FIT. 

3 

THE  ASSUMED  FORMAT  OF  YOUR  INPUT  FILE  IS:  (1X,1P402S.1S) 

ARRANGED  AS  X(I>  F1(X(I)>  F2(X< I) > ...  IN  EACH  RECORD, 

UITH  ONE  RECORD  FOR  EACH  X-VALUE  AND  B.C. 

DO  YOU  WISH  TO  CHOOSE  THE  UIDTH  PARAMETER? 

'YES' 

ENTER  A  VARIABLE  MESH  SCALE  FACTOR  <0.2DO,  3.000) 

0.442300 

THE  CHOSEN  AND  DEFAULTED  CONTROL  PARAMETERS  ARE! 

BASE  POINTS, +B.C  .  EPS.CONV  ,  IDQT  ,  ITMAX 

31  31  1 .OOOOOOOOOOOOOOOOE-02  3  0 

B.C.  DESCRIPTORS 

0  0  0  0 
BESIN  SMOOTH  INTERPOLATION  UITH  SSPLINESt 
MILD  INSL  ERRORS,  00  YOU  WANT  THE  DETAILS?  ('YES'  or  'NO  ') 

'NO' 

PROCEEDINO  TO  INTERPOLATION  AMYHAY  UITH  IDOT-  3  AND  IER-  34 
THE  FORMAT  OF  YOUR  OUTPUT  FILE  IS*.  (IX,  1P3023. 13) 

ARRANGED  A8  X(I>  ,  0(1)  ,  Ll( I)  ,  L2(I>...  IN  EACH  RECORD 
UITH  A  HEADER  FILEt  01  TITLE  BOM  DIM... ADDED 
FILES  CLOSED. 

Directory  _DRA1 JCIPR1 .NET NODE 3 

LAHZTB19 • DAT  1 1  S/ll  2A-APR-1782  14)42  (RUED. RUED, RUE, R) 

LAHZTB1PH.DAT) 1  1/11  26-APR-1982  14)42  (RUED, RUED, RUE, R> 

Total  of  2  file*,  ?/22  block*. 

■st 


Th#  two  output  file*  contain  the  ba*o  point*,  width  r*r*oot*r«,  and 
fit  coefficients  ~ 
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TYPE  COWWAMD85CHEKIT.COM 


««  EVALUATE  FILENAME  »>> 

FI 

ASSIGN  'Ft '.DAT  INFILE 
ASSIGN  'Ft 'H. DAT  HEADER 

INQUIRE  P2  *UHAT  NAME  DO  YOU  WISH  FOR  THE  OUTPUT  FILE?  ENTER  IT. PLEASE.* 
ASSIGN  'P2' .DAT  OUTFILE 
ASSIGN/USER-MODE  S YS (COMMAND i  SYStINFUT  5 
RUN  UTILITYtCHCKIT.EXE 

<««  PROMPTS  HILL  DE  8ENERATED  FOR  FURTHER  INFORMATION.  »» 

DCASSION  INFILE 
DEARS I ON  OUTFILE 
DCASSION  HEADER 
FILES  a.  DAT 


and  Its  mss  is  assin  indicated  in  the  coaaand  fils. 

EVALUATE  LAHZTB19 

WMT  NAME  DO  YOU  M1SH  FOR  THE  OUTPUT  FILE?  ENTER  IT. PLEASE. 5  BRANCH. 

THE  INTERPOLATION  FILE  FORMAT  APPEARS  TO  BE 
( IX. 1P3D23. 13) 

THE  DEFAULT  SPLINE  CHOICE  IS 
SSPL1NCS5 

ENTER  XHIN.XHAX. NUMBER  OF  INTERIOR  CHECK  POINTS: 

“3 •  3A333S7 A 1 44SSS t DO  2.333778937S44778D0  49 

DO  YOU  NISH  TO  CHANGE  THE  SPLINE  CHOICE? 

'NO' 

ENTER  THE  INDEX  OF  THE  INTERPOLANT.  NUMBERS  OF  DERIVATIVES  (0-4)5 

t  2 

ACTIVE  CONTROL  I  DOMAIN  PARAMETERS: 

XHIN  XNAX  JORD  OUTPUT  FORMAT 

-3.3433387*1440381  2 . 33377S93784477S  1  ( IX. 1P4E18.4) 

DO  YOU  UI3N  TO  CALCULATE  OTHER  INTERPOLANTS  IN  THIS  FILE? 

'NO'-* 

FILES  CLOSED. 

Director*  _DRA1 5C1PR1 .METN0DE3 

BRANCH. DAT i 1  8/11  24-APR-1982  13507  (RUED. RUED. RUE. R) 

LAMZTB19.DAT II  8/11  24-APR-1982  13507  (RUED. RUED. RUE. R) 

LAHZTB19H.DAT! 1  1/11  24-APR-1982  13507  (RUED. RUED. RUE. R) 

Total  of  3  files.  17/33  blocks. 


Here  the  output  file  contains  the  evaluation  arid,  the  inter- 
rolant.  and  its  first  two  derivatives  which  can  be  used  to 
ehecfc  the  euslitu  of  the  fit  asainst  the  orislnal  data. 

twee  BRANCH. DAT »1 
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/.v 


-3.3433396+00 

1 >0398326+00 

1.22798SE+00 

-3.239513E+00 

-3  «  4473146+00 

1.1410746+00 

8.4229026-01 

-3.1403106+00 

-3 • 3292946+00 

1.2400396+00 

3.0997836-01 

-2.4487756+00 

-3.2112726+00 

1.2030006+00 

2.470824E-01 

-1.4374416+00 

-3.0932996+00 

1.3079076+00 

1.3029396-01 

-5.5799086-01 

-2.9732276+00 

1.3233756+00 

1.4743336-01 

1.3078386+00 

-2.0372056+00 

1.3349336+00 

4.3439316-01 

3.0214426+00 

-2.7391026+00 

1.4290046+00 

7.3843846-01 

1 .4474306+00 

-2.4211006+00 

1.3103426+00 

4.4143486-01 

-4.1973436+00 

-2.3031306+00 

1.3310306+00 

-1.7983446-01 

-8 . 4472B1E+00 

-2.3031156+00 

1.4799006+00 

-8.8990826-01 

-1.8132826+00 

-2.2470936+00 

1.3019046+00 

-4.4487836-01 

3.1482226+00 

-2.1490716+00 

1.3104776+00 

-4.8349726-01 

-3.3417846+00 

-2.0310406+00 

1.2043716+00 

-1.042440E+00 

-1.043433E+00 

-1.9130246+00 

1.0837376+00 

-9.073224E-01 

2.478937E+00 

-1.7950036+00 

9.9541036-01 

-3.9309986-01 

2.338012E+00 

-1.4749016+00 

9.3774716-01 

-4.3709496-01 

-1 .3430396-01 

-1.3509396+00 

0.7748316-01 

-4.3319846-01 

-2.4302976+00 

-1.4409346+00 

7.8403936-01 

-8.8044796-01 

-1 .031771E+00 

-1.3229146+00 

4.8029436-01 

-8.7427086-01 

8.3731426-01 

-1.2040926+00 

3.8404816-01 

-7.3240346-01 

1.0437296+00 

-1.0040496+00 

3.0240496-01 

-4.311419E-01 

9.9723276 -01 

-9.4004706-01 

4.3472436-01 

-3.1924436-01 

8.9390886-01 

-0.3002446-01 

3.7947326-01 

-4.1840926-01 

8.2241086-01 

-7.3200236-01 

3.3347136-01 

-3.2527716-01 

7.490377E-01 

-4.1477996-01 

3.0222716-01 

-2.4413806-01 

4.1344116-01 

-4.9473736-01 

2.7727426-01 

-1.8227746-01 

4.34132SE-01 

-3.7073326-01 

2.3842126-01 

-1.3990936-01 

2.9740406-01 

-2.4071206-01 

2.4301196-01 

-1.0873326-01 

2.4139946-01 

-1.4249036-01 

2.3238096-01 

-8.211843E-02 

2.0423746-01 

-2.4440136-02 

2.2418746-01 

-4.1742086-02 

1.3047806-01 

9.3334226-02 

2.1734006-01 

-3.2012796-02 

1.707417E-02 

2.1137446-01 

2.1117046-01 

-S.7714406-02 

-9 . 472043E-02 

3.2939096-01 

2.0347916-01 

-7.4111146-02 

-1.7205406-01 

4.4742136-01 

1.9344346-01 

-9.3827346-02 

-1.8444076-01 

3.4344346-01 

1.8093396-01 

-1.1320286-01 

-1.349304E-01 

4.0344406-01 

1.4440106-01 

-1.2439836-01 

-5.330437E-02 

0.0140046-01 

1.3148046-01 

-1.2834936-01 

1.4989076-02 

9.1931076-01 

1.3434406-01 

-1.2419776-01 

3.0949406-02 

1.0373336+00 

1.2229426-01 

-1.1478346-01 

7.449413E-02 

1 . 1333336+00 

1.0908296-01 

-l .0475036-01 

9.3045476-02 

1.2733706+00 

9.7132806-02 

-9.390212E-02 

8.5823386-02 

1.3914006+00 

0.4344336-02 

-0.4928216-02 

4.8923716-02 

1.3094226+00 

7.4394446-02 

-7.833483E-02 

8.137033E-02 

1.4274436+00 

4.7959736-02 

-4.787143E-02 

8.8249476-02 

1.7434476+00 

4.0402206-02  ' 

-3.9738886-02 

4.3311446-02 

1.0434906+00 

3.3430816-02 

-3.4404886-02 

2.4493286-02 

1.9017126+00 

4.7234976-02 

-3.0200546-02 

9.1434486-02 

2.0997341+00 

4.2070926-02 

-3.7290246-02 

1.003034E-01 

2.2177376+00 

3.0074326-02 

-3.3841446-02 

-4.394240E-02 

2.3337796+00 

3.3193876-02 

-3.1900196-02 

-2.1410746-01 

Tha  original 

data  for  this  chack 

i*  tha  tacond 

coloan  of 

CRC19L2TB  I • E •  tha  natural  loaarith*  of  tha  branching  ratio. 


•o«  CHCMM0DCICRC19LZTB.DAT 


APPENOIX  II.  Derivation  of  the  Electrodiffusive  Limit 

The  general  solution  of  the  Hertz  vector  wave  equation  is  a  super¬ 
position  of  incoming  and  outgoing  cylindrical  waves.  The  linear  or  nonlinear 
response  of  the  medium,  specified  by  the  convolution  involving  3  J,  acts  to 

T 

cut  off  these  freely  propagating  components.  One  may  view  this  process  as 
the  production  of  reflected  power  of  nearly  opposite  intensity  to  that  incident, 
leaving  a  residual  amount  deposited  in  the  load  as  Joule  heating.  Because 
of  this  near  cancelation  of  Z  field  components  it  is  possible  to  neglect 
the  second  order  time  derivative  in  the  Z  wave  equation  in  a  first  approxima¬ 
tion.  In  particular  the  requirement  of  a  strict  detailed  balance  in  the  free 
wave  component  E*(x«t)  demands  that  72E*  *  0,  because  one  is  interested  in 

the  limit  for  which  32E*«3  J  is  true  to  arbitrary  precision.  Physically  the 

re  r  32Ew 

common  situation  is  one  of  nearly  balanced  waves  in  which  e  =  sup  «1, 

,  w  T 

but  the  free  wave  components  are  solutions  of  Q2  E2  *  0.  The  constraint 
that  72E*  *  0  is  simply  a  way  to  select  the  limit  e  -*■  0.  Since  any  solution 
of  Laplaces  equation  valid  inside  some  radius,  and  vanishing  on  that  radius, 
must  vanish  everywhere  within,  then  in  the  limit  of  truly  detailed  balance 
Inside  a  radius  xQ  and  on  that  radius,  x"*3  Cx3  E?)  *  0  is  a  necessary  condi- 

tion.  But  using  Faraday's  law  one  has  x”’*3  (x3_8a)*0  as  well.  In  terms  of  the 

X  >■  0 

Hertz  vector  Z,  8*  »  -3  3  Z,  so  that  the  original  condition  of  detailed 
8  t  x 

balance  in  the  wave  components  becomes  -x"*3,  Cx3  3iZ)  *  0  for  some  region 

X  X  1 

X  <  x0-  The  equation  governing  the  electodynamics  in  the  limit  of  strict 
detailed  balance  is  thus  obtained  by  operating  on  the  original  wave  equation 
with  72  and  removing  the  contribution  arising  from  32Z.  The  application  of 
72  on  the  convolution  with  3tJ  produces  just  the  integrand 
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\  [  V*  «roT-»  *•  £ *rx  <*»«  £  v]} 


4ir3T  JCx.t). 


The  resulting  equation  is 


X-13  x  3  E  -  4it3  J  , 
X  X  t 


(A.II.l) 


because  -V2Z  ■  E  in  the  original  relationship.  It  is  this  relationship 
which  can  be  transformed  into  a  diffusion  equation  for  E2(x>t). 

Using  the  usual  Ohm's  law,  one  finds 


V  -£1,1*11,  (Ez  ♦  Sr8e  ♦  Eth) 


E  3  I  +  1(3  E  +  S  3  8a  +  BJ  8  +3  Etu) 
T  v  T  z  Br  T  0  0  T  r  T  th; 


and  applying  Faraday's  law  once  again,  to  eliminate  3_BQ  and  form  the 

l  U 

(dimensionless)  material  derivative. 


V  -  E  »t  r  *  1  (fe  Ez  +  Wr  +  3rEth>- 


(A. II. 2) 


Here  £-  ■  3  +  S  3  ;  E  =  E,  +  8  8,  +  E...  ,  the  total  electric  field  in 
ut  t  r  x  2  r  9  tn 

r0 

the  plasma  load;  and  Z  *  a  ~  is  the  dimensionless  conductivity.  Finally, 

In  this  case  of  a  linear  response  relating  J  and  E,  the  equation  determining 
the  diffusion  of  E.  into  the  plasma  becomes 


hh  X’13JX  3  E)  -  EVnZ  - 


3t  Eth* 


(A. II. 3) 


In  the  nonlinear  case,  reflecting  a  marginally  stable  drift  speed 

limit,  this  derivation  produces  the  same  result,  insofar  as  J  a  1(E)  E  is 

formally  admissible  as  a  constitutive  relation.  However,  one  is  forced 

either  to  neglect,  or  to  model  from  first  principles,  the  term  3^£n£  in  the 

relationship  just  given.  The  reason  stems  from  the  implicit  definition  of 
en  c 

2(E)  *  cTT -  •  If  an  explicit  time  derivative  is  to  be  calculated  for 

E*Jscale 

£,  then  it  must  arise  through  some  model  of  the  turbulent  fields  providing 
the  enhanced  drag  on  one's  proverbial  test  particle.  This  model  would 
necessarily  involve  various  functionals  of  E  and  its  time  history  and  would 
generally  not  hold  £  to  the  exact,  marginally  stable  value  at  all  times. 

Such  a  model  Is  beyond  the  present  scope  of  this  work. 

It  is  a  better  choice  to  assume  that  £(E)  has  no  well-defined  time 
derivative,  that  it  respresents  a  local,  time  averaged  estimate  of  the  truly 
physical  transport  coefficient.  A  strict  application  of  the  marginal 
stability  constraint  decouples  J  and  E  once  E  becomes  sufficiently  large. 

When  this  decoupling  is  Imposed  on  a  diffusion  equation  there  is  no  mechanism 
left  to  evolve  J  and  E,  cf.  Eqn(A.II.l).  The  more  accurate  first  principles 
model  would  not  force  a  complete  decoupling,  but  would  provide  a  smooth 
transition  between  linear  response  and  drift  speed  limitation  —  the 
appropriate  asymptotic  conditions.  In  contrast,  for  the  complete  electro¬ 
magnetic  description,  a  drift  speed  limit  is  fully  admissible  —  it  simply 
specifies  a  wave  source  term  »  e3T(necs)*  The  use  of  the  full  wave 
equation  is  preferable  whenever  strict  adherence  to  the  drift  speed  limitation 
is  desired,  and  the  wave  components  generated  in  the  drift  speed  limited 
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regions  of  the  load  are  probably  important  contributions  to  the  dynamics  of 
the  transition  front  between  the  two  asymptotic  domains  of  resistivity.  On 
the  other  hand.  If  the  diffusive  approximation  is  adopted,  then  the  drift 
speed  limitation  must  be  applied  with  the  caveat  that  its  time  dependence 
3TinZ(£)  cannot  be  specified  from  the  local  drift  speed  constraint  alone 
and  is  better  neglected  than  modeled  in  any  ad  hoc  fas ion. 


APPENDIX  III.  Area  Weighted  Differences 

In  the  Lagranglan  mesh  shown  in  Figure  4  the  only  fundamental 

positions  are  the  radial  cell  boundary  locations  *ry  and  their  associated 

velocities  *r..  The  scalar  quantities  assigned  to  the  cell  centers  have 
J 

Implied  spatial  locations,  *r.,  derived  from  the  boundary  locations.  In  order 

C  J 

to  avoid  formal  singularities  in  various  derivative  expressions  the  independent 
variable  is  chosen  here  as  a^  ■  C**j)  »  an  areal  coordinate.  The  two 
primary  types  of  differential  operators  to  be  estimated  are  CD  derivatives 
of  first  and  second  order  at  a  boundary  based  on  cell-cenetered  or  cell 
boundary  data,  and  (il)  derivatives  at  the  cell  centers  based  on  cell- 
centered  data. 

Turning  first  to  those  operators  based  on  or  producing  cell  boundary 
data,  the  radial  gradient  for  any  function  F.  defined  on  *r.  is  estimated  by 

C  J  C  J 


O  F)  .  *  2  *r  , 
r  J  J\caj  cVl 


(A. Ill .1) 


with  the  second  derivative  estimate,  based  on  F.,  given  by 

J 


Of  F),*2«y,  2*r 


F.^,  -  F.' 


-  2*r , 


r  s  jV3Wr»ii  c  ^\ai 

Caj  '  caj-l 


,  (A. III. 2) 


and  the  averaging  operator  to  the  boundary  defined  as 


a.  -  a 


F.  *  s-* - - — 

J  \caj  *  caj-l 


p  )  F 

c  J"1  \caj  -  caj-l/  c  j 


(A. III. 3) 


These  difference  operators  are  then  sufficient  to  cover  all  gradients  and 
averages  needed  in  computing  the  evolution  of  the  fundamental  boundary 


variables  [Vr ( t ) ] ^  «  *r^  and  Its  integral,  the  Lagrangian  position. 
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The  corresponding  operators  estimating  first  and  second  order 


gradients  at  the  cell  centers  can  be  obtained  similarly.  Here  the  use  of 

ghost  points  located  at  r  and  r.^,  is  the  natural  means  of  enforcing  the 

c  o  c  J+l 

boundary  conditions  appropriate  to  the  diffusive  processes,  such  as  thermal 
condition.  The  difference  operators  given  below  assume  that  some  cell -centered 
quantity  cFj  has  been  assigned  ghost  point  values  cFj+j  *  c ,  cFq  *  CF^  in 
order  to  make  the  flux  estimate,  proportional  to  3rF,  vanish  at  the  first  and 
last  boundary.  First  one  may  define  several  Cnonuniform)  mesh  intervals. 


hj  ■  caj  -  caj-l 


hj  ■  cVl  -  caj 


CA.m.4) 


j  j, 

and  the  ghost  point  locations 


h0>*  2<cal  •  al> 


2(aJ+l  -  Ca0>  • 


The  first  derivative  3  F  is  then  estimated  by 


c(3aF)l* 


a 

rh< 

JL 

F i  *  h^  F.  -  -J-  F,  , 

•  * 

hi 

C  j+l  j  Crj  h<  c  j- 1 

CA .111.5)  £ 

with  boundary  conditions. 


,  p)  ho  (cFZ~cFV 


'  C<V>N 


fcFj'cFJ-l 


hJ+hJ-l 


J-l 


The  second  derivatives  are  similarly  obtained  using 


S&s  *  ;*v 

hj  hj 


"j  hj 

*7  ^  ^  +  *3  cFj-l 


while  the  boundary  conditions  become 


2_ 

> 


cFj"cFj-r 


'j-i 


^-i+hj 


(A. III. 6) 


All  the  expressions  reduce  to  well  known  second  order  accurate 


differencing  schemes  when  the  mesh  is  presumed  uniform. 
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Lawrence  Livermore  Lab 
ATTN:  FC-1 

Field  Command 

Defense  Nuclear  Agency 
ATTN:  FCPR 
ATTN:  FCTT,  W.  Sumna 
ATTN:  FCTT 
ATTN:  FCTXE 
ATTN:  FCTT,  G.  Ganong 


Under  Secy  of  Def  for  Rsch  A  Engrg 

ATTN:  Strategic  A  Space  Sys,  OS 


DEPARTMENT  OF  THE  ARMY 


Harry  Olanond  Laboratories 
ATTN:  DELHD-NW-P 
ATTN:  DELHO-NW-RA 
ATTN:  DELHD-NW-RI,  Kervls 
ATTN:  OELHO-TA-L 


DEPARTMENT  OF  THE  NAVY  (Continued) 

Naval  Surface  Weapons  Center 
ATTN:  Code  R40 
ATTN:  Code  F31 
ATTN:  Code  F34 

DEPARTMENT  OF  THE  AIR  FORCE 

Air  Force  Weapons  Laboratory 
ATTN:  SUL 
ATTN:  CA 
ATTN:  NT7C 
ATTN:  NT 
ATTN:  HTTP 

Ballistic  Missile  Office 
ATTN:  ENSN 
ATTN:  SYDT 

Deputy  Chief  of  Staff 

Research,  Development,  A  Acq 
ATTN:  AFROQI 

Space  Division 

ATTN:  XR,  Plans 
ATTN:  YEZ 
ATTN:  YGJ 
ATTN:  YKF 

ATTN:  YKS,  P.  Stadler 
ATTN:  YKM 
ATTN:  YNV 

Strategic  Air  Comand 

ATTN:  I NT,  E.  Jacobsen 

DEPARTMENT  OF  ENERGY 


Office  of  Military  Application 
Attention  Ofc  of  Inert  Fusion  for 
ATTN:  S.  Kahalas 
ATTN:  T.  Godlove 
ATTN:  C.  Hllland 


US  Army  Nuclear  A  Chemical  Agency 
ATTN:  Library 

US  Army  Test  and  Evaluation  Comd 
ATTN:  DRSTE-a 

USA  Missile  Coanmnd 

ATTN:  Oocuaents  Section 


Naval  Research  Laboratory 

ATTN:  Code  2000,  J.  Brown 

ATTN:  Code  4770,  I.  Vltokovltsky 

ATTN:  Code  4720,  J.  Oavls 

ATTN:  Code  4780,  S.  Ossakow 

ATTN:  Code  4773,  G.  Coopers  tel  n 


Naval  Weapons  Center 

ATTN:  Code  343,  FKA6A2,  Tech  Svcs 


OTHER  GOVERNMENT  AGENCY 

Central  Intelligence  Agency 
ATTN:  OSWR/NED 


DEPARTMENT  OF  ENERGY  CONTRACTORS 


University  of  California 
Lawrence  Livermore  National  Lab 


ATTN: 

ATTN: 

ATTN: 

ATTN: 


Technical  Info  Dept  Library 
L-545,  J.  Nuckolls,  Class  L-33 
L-47,  L.  Wouters 
L-153 


ATTN: 

ATTN: 


L-153,  D.  Meeker,  Class  L-477 
W.  Pickles,  L-401 


Los  Alamos  National  Laboratory 
ATTN:  MS222,  J.  Brownell 


DEPARTMENT  OF  ENERGY  CONTRACTORS  (Continued)  DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued! 


Sandla  National  Lab 

ATTN:  M.  Clauser,  Org  5241 
ATTN:  6.  Kuswa,  Org  5240 
ATTN:  S.  Yonas 
ATTN:  Tech  Lib  3141 
ATTN:  0.  Powell 
ATTN:  Org  9336,  0.  Allen 

DEPARTMENT  OF  OEFENSE  CONTRACTORS 


Advanced  Research  &  Applications  Corp 
ATTN:  R.  Artni stead 


Aerospace  Corp 

ATTN:  V.  Josephson 
ATTN:  S.  Bower 

ATTN:  Technical  Information  Services 
BOM  Corp 

ATTN:  Corporate  Library 

BOM  Corp 

ATTN:  L.  Hoeft 

Boeing  Co 

ATTN:  Aerospace  Library 
Dlkewood 

ATTN:  Tech  Lib  for  L.  Davis 

EGAG  Nash  Analytical  Svcs  Ctr,  Inc 
ATTN:  Library 


General  Electric  Co 
ATTN:  J.  Peden 
ATTN:  H.  O'Donnell 

IRT  Corp 

ATTN:  R.  Hertz 

JAYCOR 

ATTN:  E.  Wenaas 
JAYCOR 

ATTN:  E.  Alcaraz 
ATTN:  R.  Sullivan 
4  cy  ATTN:  R.  Terry 
4  cy  ATTN:  J.  Guillory 

Kaman  Sciences  Corp 
ATTN:  S.  Face 


Lockheed  Missiles  &  Space  Co,  Inc 
ATTN:  L.  Chase 

Lockheed  Missiles  &  Space  Co,  Inc 

ATTN:  S.  Talmuty,  Dept  81-74/154 

Maxwell  Labs,  Inc 

ATTN:  A.  Kolb 
ATTN:  D.  Tanlmoto 
ATTN:  A.  Miller 
ATTN:  0.  Cole 

McDonnell  Douglas  Corp 

ATTN:  S.  Schneider 

Mission  Research  Corp 
ATTN:  U.  Hart 
ATTN:  C.  Longmlre 

Mission  Research  Corp 
ATTN:  B.  Godfrey 

Mission  Research  Corp,  San  Diego 
ATTN:  B.  Passenhelm 

Pacific-Sierra  Research  Corp 

ATTN:  H.  Erode,  Chairman  SAGE 
ATTN:  L.  Schlesslnger 

Physics  International  Co 
ATTN:  C.  Gilman 
ATTN:  C.  Stallings 
ATTN:  G.  Frazier 

RAD  Associates 

ATTN:  P.  Haas 
ATTN:  A.  Latter 

RAD  Associates 

ATTN:  P.  Turchl 

S-CUBED 

ATTN:  A.  Wilson 

Science  Applications,  Inc 
ATTN:  K.  Sites 

TRW  Electronics  A  Defense  Sector 

ATTN:  Technical  Information  Center 
ATTN:  D.  Clement 


Kaman  Tempo 

ATTN:  DAS I AC 


END 
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